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Abstract 

Single  surface  multipactor  is  a  phenomenon  which  imposes  a  limit  on  the  total  power 
that  may  be  transmitted  by  a  High  Powered  Microwave  (HPM)  System.  Electron 
multipactor  results  in  an  avalanche  of  high  energy  electrons  within  the  waveguide.  In 
a  relatively  short  time,  the  electrons  deposit  power  on  the  surface  of  the  dielectric 
window  which  separates  the  waveguide  from  the  antenna  and,  eventually,  the  window 
can  be  damaged  and  the  system  will  fail.  Mitigation  approaches  for  single  surface 
multipactor  at  dielectric  windows  are  investigated  using  Particle-In-Ccll  (PIC)  sim¬ 
ulations.  Initially  baseline  susceptibility  diagrams  are  constructed  analytically  and 
compared  with  self-consistent,  dynamic  system  trajectories.  The  power  deposited 
on  the  surface  of  a  dielectric  window  in  an  HPM  system  is  considered  using  three 
different  methods  and  the  results  of  PIC  simulations.  Geometric  mitigation  is  then 
considered  by  varying  the  window  orientation  with  respect  to  the  HPM  electric  field. 
Small  angular  deviations,  less  than  20  degrees,  from  the  nominal  case  of  normal  in¬ 
cidence  show  dramatic  changes  in  the  susceptibility  diagram.  A  materials  approach 
to  mitigation  is  then  considered.  Titanium  Nitride,  TiN,  coatings  applied  to  the 
dielectric  surface  can  substantially  reduce  the  secondary  emission  yield.  Represen¬ 
tative  modifications  of  the  secondary  emission  yield  are  simulated  and  the  resulting 
susceptibility  diagrams  are  discussed. 
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ELECTRON  MULTIPACTOR:  THEORY  REVIEW,  COMPARISON 
AND  MODELING  OF  MITIGATION  TECHNIQUES  IN  ICEPIC 


I.  Introduction 


1.1  Problem  Description 

Recent  technological  advances  in  the  study,  design  and  implementation  of  High 
Powered  Microwave  (HPM)  Devices  have  led  to  increased  power  levels  throughout  the 
entire  system.  Higher  power  levels  present  many  problems  within  the  waveguide  that 
transfers  the  energy  between  the  source  and  the  antenna,  most  notably  breakdown 
and  electron  multipactor.  In  order  to  prevent  breakdown,  the  waveguide  must  be 
kept  at  very  low,  near-vacuum  pressures.  A  dielectric  window  situated  at  the  end 
of  the  waveguide  maintains  the  low  pressure  inside  the  waveguide,  while  allowing 
the  Radio  Frequency  (RF)  energy  to  pass  through  to  the  antenna.  Single  surface 
multipactor  can  occur  at  the  surface  of  this  window;  which,  when  combined  with 
thermal  heating  due  to  the  high  power  throughput,  can  lead  to  cracks  in  the  window 
and,  soon  thereafter,  catastrophic  system  failure.  Therefore,  it  is  desirable  and  even 
necessary,  for  future  systems,  to  mitigate  this  phenomenon. 

1.2  Background 

Electron  multipactor  is  a  phenomenon  first  observed  in  1934  by  Otto  Farnsworth 
[7].  Farnsworth  was  able  to  harness  the  effect  to  create  an  ”AC  electron  multiplier,” 
thereby  terming  the  phenomenon  ”  multipactor” .  The  following  brief  description  of 
electron  multipactor  will  be  followed  by  a  more  detailed  discussion  in  Chapter  II. 
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Multipactor  occurs  in  two  distinct  cases:  in  a  waveguide  or  RF  cavity  and  at  a  single 
dielectric  surface.  In  both  cases,  ’’seed”  electrons  are  provided  by  stray  electrons 
emitted  by  random  emission  processes  (cosmic  rays,  material  emission,  etc).  In  the 
first  case,  microwave  energy  propagates  through  a  (semi-infinite)  waveguide  cavity, 
an  AC  electric  field  is  created  between  the  two  sides  of  the  cavity  and  electrons  from 
one  side  are  accelerated  towards  the  other  side.  When  these  energetic  electrons  strike 
the  surface,  more  than  one  electron  may  be  emitted  according  to  Vaughn’s  formula 
for  secondary  emission  [17].  This  leads  to  a  loading  and  detuning  of  the  cavity.  In 
the  second  case,  at  the  surface  of  a  single  dielectric,  electrons  leave  the  surface  of 
the  dielectric,  leaving  behind  an  overall  positive  surface  charge.  These  electrons  then 
gain  energy  from  the  microwave  field,  but  are  returned  to  the  surface  by  the  restoring 
field  generated  by  the  separated  charge  in  the  system  (referred  to  as  EDC ).  As  in 
the  previous  case,  upon  striking  the  dielectric  surface,  more  than  one  electron  may  be 
emitted.  Overall,  as  more  and  more  electrons  are  emitted  and  strike  the  surface,  they 
will  deposit  a  certain  amount  of  power  on  the  surface,  which  can  lead  to:  destruction  of 
the  dielectric  material  and,  in  the  case  of  High  Powered  Microwave  (HPM)  systems, 
loss  of  necessary  vacuum  conditions,  perpetuation  of  further  loss  mechanisms  and 
eventual  failure  of  the  system.  According  to  Barker  [3],  multipactor  is  not  the  sole 
cause  of  window  failure,  rather,  it  is  an  enabling  process  in  a  series  of  destructive 
phenomena  that  occur  under  some  circumstances  in  a  single-surface  geometry. 

While  general  theory  attracted  a  small  amount  of  attention,  not  until  the  advent 
and  widespread  use  of  modern  microwave  devices  (especially  those  at  higher  powers) 
was  a  coherent  theory  formulated  by  Vaughn.  The  culmination  of  his  work  [17]  is  the 
foundation  of  multipactor  theory  and  is  commonly  referenced  to  describe  secondary 
emission  from  any  given  material  using  a  small  number  of  parameters.  In  the  past 
decade,  a  small  group  of  researchers  (Ang,  Lau,  Verboncoeur,  Kishek,  Sanzatov, 
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Valfclls,  etc.)  have  continued  to  investigate  and  develop  theoretical  representations 
of  multipactor.  Their  work  is  summarized  in  Chapter  II. 

1.3  Mitigation  Concepts 

Mitigation  concepts  have  been  investigated  in  the  past,  but  conclusive  data  is  still 
lacking  in  the  effectiveness  of  such  strategies.  Two  strategies  are  considered  in  this 
paper: 

•  Thin  Film  Coatings  -  a  thin  coating  of  a  material  with  a  low  Secondary  Electron 
Yield  (SEY)  ,  ideally  at  or  near  one,  is  deposited  on  the  surface  of  the  dielectric, 
with  the  goal  of  reducing  the  amount  of  electrons  produced  by  multipactor  in 
the  system. 

•  Geometrical  Modification  of  the  Waveguide/Window  Interface  -  an  angle  of 
obliqueness  is  introduced  between  the  orientation  of  the  RF  held  (Erf)  and  the 
dielectric  surface,  resulting  in  a  narrowing  of  the  susceptibility  of  the  system  to 
multipactor. 

1.4  Research  Goals 

The  goal  of  this  thesis  is  twofold:  to  investigate  previous  research  efforts,  validating 
both  analytical  and  computational  results  of  those  efforts  and  to  seek  out  methods  to 
prevent  electron  multipactor  or  mitigate  the  effects  thereof,  in  the  context  of  HPM 
systems.  The  two  primary  methods  considered  for  mitigation  are:  thin  coatings 
of  special  materials  and  geometric  modifications  to  the  window.  Thin  coatings  of 
materials  with  lower  secondary  emission  coefficients,  Smax o,  reduce  the  overall  quantity 
of  secondary  electrons.  Under  some  circumstances,  this  method  may  prevent  the 
avalanche  of  secondaries  and  the  deposition  of  power  that  results  in  damage  to  the 
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dielectric  window.  The  properties  of  these  materials  will  be  further  discussed  in 
Chapter  II.  The  geometric  modification  method  relies  on  the  interaction  of  the  electric 
fields  near  the  surface  of  the  dielectric  surface  to  alter  the  susceptibility  of  a  system 
to  multipactor,  a  result  that  will  also  be  discussed  in  detail  in  Chapter  II. 
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II.  Multipactor  Theory 


In  order  to  further  understand  multipactor  and  to  analyze  possible  mitigation 
strategies,  the  analytical  theory  must  be  well  understood.  Current  theory  on  multi¬ 
pactor  is  well  described  and  documented  by  a  small  group  of  researchers,  including, 
but  not  limited  to:  Verboncoeur,  Kishek,  Lau,  Ang,  Fichtl  and  Kim.  A  summary  of 
their  work  and  conclusions  is  presented  in  this  chapter. 

2.1  RF  Cavity  Multipactor 

While  a  lengthy  discussion  of  multipactor  in  a  cavity  would  be  productive  and  in¬ 
formative,  the  primary  focus  of  this  research  is  on  single  surface  multipactor.  However, 
for  completeness,  a  short  discussion  of  the  theory  associated  with  cavity  multipactor 
follows.  The  effect  of  multipactor  in  a  cavity  system  is  well- described  by  Kishek  and 
Lau  [10].  They  compare  cavity  multipactor  to  an  equivalent  RF  circuit,  in  which  cav¬ 
ity  loading  and  detuning  are  the  primary  loss  mechanisms.  The  generic  circuit  they 
assume  is  given  in  Figure  1,  where  a  single  sheet  of  electrons  of  density  a  traverses  the 
cavity  of  width  D  as  time  progresses,  leading  to  a  (normalized)  multipactor  current 
of: 


T  .  ,  .  dx(t) 

Im  (t)  =  -a  (t) 

This  current  is  assumed  to  be  the  current  which  loads  the  cavity  and,  therefore, 
diminishing  the  RF  field  as  it  propagates  through  the  system.  This  acts  as  a  loss 
mechanism,  reducing  the  total  power  output  of  the  HPM  system. 

The  secondary  emission  of  the  two  surfaces  is  exactly  the  same  as  that  for  the 
dielectric  surface  (given  in  further  detail  in  Section  2.2).  It  is  determined  by  Kishek 
and  Lau  [10]  that  the  RF  beam  loading  effects  are  much  more  important  than  the 


5 


Id  I 


1 1  1 

+  Inv 

D 

F)  i  <  £3 

L  ]vg  U  V(0t 

Y  '•  F3 

L  ' 

c:  L.. 

i  ?  F 

■ .  *(«> 

RF  Cavily 


Mullipactor 


Figure  1.  Kishek  and  Lau’s  equivalent  RF  circuit  [10] 


space  charge  effects  for  a  high  quality  factor,  Q  =  ^  y  §  >  10,  and  for  a  small  spatial 
extent,  x,  the  power  delivered  to  the  surfaces  by  multipactor  can  be  very  high.  This  is 
a  simplistic  model,  designed  only  to  show  that  multipactor  can  indeed  deliver  power 
to  the  walls  of  the  cavity  and  reduce  the  output  power  of  the  system. 


2.2  Single  Surface  Multipactor 

Single  surface  multipactor,  such  as  that  occurring  at  the  surface  of  the  dielectric 
window  in  an  HPM  device,  is  a  mechanism  resulting  in  power  deposition  on  the 
surface  of  the  dielectric,  physical  damage  to  the  dielectric  window  and  catastrophic 
failure  of  the  HPM  system.  The  phsyical  consequence  of  this  effect  is  a  limit  in  the 
amount  of  power  that  may  be  transported  through  the  waveguide  and  to  the  target. 
It  is  also  worth  noting  that  electron  multipactor  acts  as  a  beam-loading  mechanism, 
decreasing  the  total  power  output  of  the  system. 

The  geometry  for  single  surface  multipactor  is  given  by  Figure  2.  Since  the  goal 
of  this  research  is  to  reduce  the  amount  of  power  deposited  on  the  dielectric  surface, 
it  is  worth  noting  that  early  work  by  Kishek,  et  al.  [10]  estimated  the  power  de¬ 
posited  on  the  dielectric  window  to  be  on  the  order  of  1%  of  the  total  RF  energy. 
This  estimate  will  be  compared  with  the  later  analytical  developments  of  multipactor 
theory.  Kishek,  et  al.  were  the  first  to  generate  susceptibility  curves  for  a  wide 
range  of  materials.  Early  developments  of  multipactor  theory  relied  on  a  simple,  yet 
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unrealistic  assumption  that  neglected  space  charge  effects.  For  the  purposes  of  this 
research,  Space  Charge  (SC)  effects  refer  to  the  electromagnetic  contributions  of  elec¬ 
trons  to  the  system.  Therefore,  when  space  charge  is  ignored  (called  the  Non-Space 
Charge  (NSC)  case),  the  current  density,  J  in  Ampere’s  Law,  is  hard-wired  to  zero. 
Physically,  this  enforces  a  stationary  solution  of  the  charge  conservation  equation  and 
neglects  the  influence  of  the  changing  electron  density.  Therefore,  the  restoring  held, 
Edc  ,  is  not  allowed  to  evolve  ”  naturally” ,  according  to  the  number  of  electrons  in 
the  system,  but  is  assumed  to  be  a  constant  value.  For  a  complete  development  of 
early  theory,  please  see  Barker  [3]  or  Fichtl  [7].  In  order  to  gain  a  more  fundamental 
understanding  of  the  complexity  of  the  multipactor  problem,  this  development  will 
start  with  a  description  of  secondary  emission,  then  move  to  an  analysis  of  the  equa¬ 
tions  of  motion  for  multipactoring  electrons  and  conclude  with  an  estimation  of  the 
power  deposited  on  the  surface  of  the  dielectric.  This  development  includes  space 
charge  effects. 

2.2.1  Conventions,  Notations  and  Definitions. 

Since  the  notational  conventions  vary  throughout  the  literature,  the  following 
standards  are  implemented  in  describing  this  work: 

•  A  block  E  represents  an  electric  held 

•  A  script  £  represents  an  energy,  such  as  the  impact  energy  of  electrons,  £% 

•  A  word  (or  words)  set  in  typewriter  font,  such  as  ENERGY,  represents  an  input 
parameter  or  section  of  input  parameters,  usually  to  simulation. 

•  Standard  physical  quantities,  such  as  mass,  acceleration  and  velocity,  are  de¬ 
noted  by  their  usual  forms,  m,  a  and  v. 
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Figure  2.  The  general  geometry  for  single  surface  multipactor 

•  Subscripts  are  used  to  distinguish  similar  quantities  in  a  system,  or  to  provide 
further  information  about  a  specific  quantity.  In  this  case,  ay  would  be  a  y- 
directed  acceleration  and  £t  would  be  the  impact  energy  of  an  electron. 

In  order  to  use  clear  and  consistent  notation  throughout  this  paper,  definitions  for 
the  terms  in  many  of  the  equations  are  as  follows  (exceptions  are  noted  and  defined 
as  they  are  encountered): 
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&max  0 
£max  (0 


=  the  ratio  of  emitted  secondaries  to  incident  electrons  as 
a  function  of  the  impact  energy  and  angle 
=  shorthand  for  5  (Si,  () 

=  the  characteristic  secondary  emission  term 
=  the  impact  energy  as  a  function  of  impact  angle 
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Smax o  =  the  maximum  impact  energy  that  produces  5max 0 

(  =  the  angle  of  impact,  measured  from  the  normal 

Si  =  the  impact  energy  of  impinging  electrons 

ks  =  a  surface  smoothness  parameter,  usually  equal  to  1 

Up  =  plasma  frequency 

u  =  RF  angular  frequency  of  the  system 


2.2.2  Secondary  Emission. 

The  most  commonly  accepted  form  for  secondary  emission  at  the  dielectric  surface 
is  given  by  Vaughn’s  model  [17]: 


*(£,0 


k 


s(Si,  0 


3max  (0 

£max  (0 


w 


Smax  (C)  '  (»e(1  w^k(w) 

{0.56  for  w  <  1 
0.25  for  1  <  w  <  3.6 
8max  (C)  •  (l.l25uT0-35)  for  w  >  3.6 


5 max 0  •  (  1  + 

Smax  0 

Sr 


1  + 


KC 

2tt 

hC 

7 T 


£ 

'-'max 


(0 


(1) 

(2) 

(3) 

(4) 

(5) 

(6) 


Where  Vaughan  later  revised  his  equations  by  adding  Equation  3  in  [18]  to  account  for 
behavior  at  large  values  of  w.  Such  large  values  of  w  are  unrealistic,  however,  in  the 
course  of  numerical  simulations,  a  few  particles  may  reach  such  values  and  must  be 
accounted  for  to  assure  stability  of  the  code.  Both  simulations  and  calculations  show 
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that  most  impact  energies  to  fall  primarily  below  the  3.6 £max  threshold.  This  is  best 
illustrated  by  Figure  3,  where  only  a  few  particles  rise  above  the  3.6£maxo  line.  Only 
when  very  small  (~  5°)  angles  of  obliqueness  are  introduced  between  Erf  and  the 
dielectric  surface  do  the  impact  energies  move  into  higher  values  of  w,  a  phenomenon 
which  will  be  discussed  in  detail  in  Section  2. 2. 5. 2. 

Valfclls,  et  al.  [16]  use  the  following  continuous  representation  for  k: 

k  =  0.435  —  - — -  arctan  ( n  In  (  —  J  J  (7) 

VT  V  \£max{QJJ 

In  order  to  better  understand  and  compare  Vaughn’s  early  and  revised  secondary 
equations  with  the  continuous  equation  given  by  Valfclls,  et  al.  [16],  Figure  4  plots 
the  three  different  forms  of  the  secondary  emission  yield  for  a  material  with  a  Smax, o 
of  3.  This  type  of  plot  allows  for  an  easier  interpretation  of  the  secondary  emission 
function.  £\  and  £2  are  commonly  referred  to  as  ’’crossover  energies”  and  represent 
the  impact  energy  at  which  the  SEY  is  exactly  one.  When  £\  <  £i  <  £2,  electron 
growth  occurs  within  the  system.  As  is  evident  from  the  piecewise  definition  of  the 
equations,  the  plot  may  be  broken  up  into  three  regions.  The  Erst  region,  which  is 
shown  in  greater  detail  in  the  inset  next  to  the  main  plot,  describes  the  SEY  of  the 
material  at  low  impact  energies  (less  than  £maxo)-  In  this  region,  the  first  crossover 
energy,  £\,  is  seen  to  differ  by  a  relatively  large  amount  when  comparing  Vaughan’s 
early  equations  with  the  continous  form  and  Vaughan’s  revised  equations.  In  the 
middle  region,  between  £maxo  (~  750  eV)  and  the  vertical  dashed  red  line  (which 
corresponds  to  an  approximate  value  of  3.6£maxo),  all  the  forms  are  seen  to  agree 
fairly  well.  In  the  region  to  the  right  of  the  ~  3.6£maxo  line,  the  forms  (especially 
Vaughan’s  revised  form)  diverge  rapidly. 
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Figure  3.  The  mean  impact  energy  of  collected  particles  at  discrete  time  steps  is  shown. 
In  this  system,  space  charge  effects  are  considered  and  the  following  parameters  are 
used:  Erf  —  5.5  MV/m,  Smax o  =  3,  £max o  =  420  eV,  /  =  2.45  x  109  GHz  and  ka  =  1.  The 
dotted  red  line  shows  the  impact  energy  which  is  equal  to  3.6£maxo-  From  this  figure,  it 
is  evident  that  very  few  particles  obtain  an  energy  requiring  Equation  3,  however,  for 
these  few  cases,  this  additional  region  is  imperative  for  numerical  stability. 


2.2.3  Kinematic  Theory. 

Using  the  relations  of  the  impact  energy  and  secondary  emission  (Equations  1  - 
7),  coupled  with  the  equations  of  motion,  the  total  energy  gained  by  the  electrons  in 
flight  may  be  determined,  along  with  their  flight  time  (also  known  as  the  ”  hoptime” 
of  an  electron)  and  the  total  power  deposited  on  the  dielectric  by  the  electrons.  The 
development  and  equations  in  this  section  are  primarily  from  Valfells,  et  al.  [16], 
where  the  equations  have  been  slightly  modified  to  match  the  conventions  described 
in  Section  2.2.1. 

Beginning  with  the  equations  of  motion,  we  can  write  the  y-directed  acceleration 
(ay)  due  to  ERF  of  an  electron  leaving  the  dielectric  (such  as  one  in  Figure  5): 


ay  — 


d2y 

dt2 


zErf 

m 


sin(cut  +  9) 


Fy 

m 


where  9  is  the  phase  of  the  RF  held  at  which  the  electron  is  emitted,  Fy  is  the  y- 
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Figure  4.  Comparison  of  the  forms  of  the  secondary  emission  equation.  The  horizontal 
dotted  line  shows  the  crossover  point  where  <5=1,  while  the  vertical  dotted  line  shows 
the  point  at  which  £t  ss  3.6£maxo. 


directed  force  on  the  electron  and  e  is  the  charge  of  the  electron.  Integrating  once 
from  0  to  r,  where  r  is  the  hoptime  of  an  electron,  we  obtain  the  y-directed  velocity: 


dy 


'Jy~  dt~  V°y  + 


eE 


RF 


mu 


(cos(ut  —  6)  —  cos{6 )) 


(8) 


where  Vqv  is  the  ^-component  of  the  electron’s  emission  velocity.  Therfore,  the  energy 
gained  by  the  electron  during  it’s  hoptime,  A£y,  can  be  determined  from 

A£>^™A vZ  =  12m((^r)  (9) 


Substituting  8  into  9,  we  find: 


A£y 

A£y 


cErf 

mu 


( cos{ur  —  6) 


cos{6 )) 


J0y 


- — [cos(ur  +  9)  -  cos(6)f 
Zmuz 

+  — ^  [, cos{ut  +  6)  —  cos{0)\ 

UJ 


(10) 

(ii) 
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X 

Figure  5.  Figure  2  is  shown  again  for  convenience.  The  general  geometry  for  single 
surface  multipactor. 


This  impact  energy  is  dependent  on  the  phase,  6 ,  of  Erf  at  which  the  particle  is  emit¬ 
ted.  Assuming  a  uniform  emission  distribution  from  0  to  2tt,  the  following  expression 
is  found  for  the  phase-averaged  energy  of  a  particle: 


(Af„}e  =  ^ 


2tt 


A  Sy  dd 


2v r 


e2E2BF  wr  e2EFF 

- sm( — )  = - ~rT  (1  —  cos(ujt)) 

men2  v  2  ’  2 mu2  v  v  ” 


(12) 


This  leaves  the  average  energy  as  a  function  of  r,  the  hoptime.  In  order  to  obtain  the 
total  average  energy,  we  now  must  average  over  the  range  of  possible  flight  times: 

r°°  p2jx;2 

(A£y)  =  (£0)  +  ^  (1  -  cos^t))  g  (r)  dr  (13) 

Where  g  (r)  is  the  distribution  function  of  all  possible  flight  times.  The  function 
g  (r)  is  dependent  on  the  distribution  function  of  electron  emission  velocities,  which 
is  usually  taken  to  be  Maxwellian.  (So)  is  the  average  emission  energy,  which  is  also 
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dependent  on  the  distribution  function  of  electron  emission  velocities.  Because  space 
charge  effects  are  included  in  this  development,  a  solution  to  Poisson’s  equation  is 
required: 


<920  (x)  en0  (  ecj)  ( x ) 

Where  0  is  the  potential  of  an  electron  density  en0,  n0  is  the  electron  density  at  x 
and  6q  is  the  free-space  permittivity.  The  number  density  profile  is  given  by: 


(14) 
=  0 


n  (x) 


vof  Qo)  dv o 
v  ( x )2 


where  v0  is  the  injection  velocity  and: 


Vmin 

v  ( x )2 

/  M 


2 ecj)  (x) 
m 


v0  +  vn 


2  n0 

V^vt 


exp 


Using  the  Maxwellian  distribution,  the  case-specific  solution  to  14  is: 


0  (x)  =  ln 

UJp 


(15) 


where  E$  is  the  surface  electric  held  and  lup  is  the  plasma  frequency.  These  are  given 
by: 


E0  = 


2kTrio 

eo 


and 


Which  yields  a  number  density,  n(x),  prohle  of: 


(  \  n° 
n  (x)  =  —  exp 

vt 


f  2 60  (g)  \ 
V  mvt  ) 


LOr) 


n0  1  4 - x 

vt 


-2 


(16) 
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Using  the  law  of  energy  conservation,  we  can  generate  an  equation  describing  the  the 
motion  of  an  electron  in  a  potential  given  by  15: 


vx  W  =  fjr  =  J' vo  ~  2vt  In  ^1  +  (17) 

By  setting  the  right  hand  side  of  this  equation  equal  to  zero,  we  can  determine  the 
maximum  excursion  of  an  electron,  since  x  =  0  describes  the  points  at  which  the 
electron  is  at  the  dielectric  surface: 


0  =>  x 


max 


(18) 


At  this  point  we  can  directly  integrate  17  to  determine  the  flight  time  for  a  given 
electron  excursion  in  the  potential  of  15: 


and,  substituting  in  18,  we  can  determine  the  maximum  flight  time  of  an  electron 
leaving  the  surface  and  returning  to  the  surface  after  making  an  excursion  of  twice 
Xmax-  This  is  the  hoptime  referred  to  in  8-13. 


r 


exp 

(u0)  =  - 


Returning  to  13,  we  can  write: 

i  roo  2  rp2 

(A£y)  =  (eQ)  +  —  RF2  (1  -  cos (u;r  (v0)))  f  M  dr 

n0  Jo  Zrnu2 
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or: 


(ASy)  —  (Eo)  + 


'•°o  ^ 2  77'2 


e  Erf  q  l  u 


I  o  2?m;2  \cnp 


Where  y  has  been  substituted  for  %,  leading  to: 


u; 


G  —  =1 


cur 


2\^2 


7T 


cos  f  exp  (_ V)  dy  (19) 


V 


y 


This  describes  the  shielding  effects  in  the  system  that  arise  from  considering  space 
charge  effects.  Now,  we  can  determine  the  injection  current  density,  Jo,  for  the 
Maxellian  distribution: 


Jo  =  e  vf  ( v )  dv  =  'f*n£eit  (20) 

From  the  impact  energy  and  the  current  density,  the  power  deposited  on  the  dielectric 
is  given  by: 


{Si)  Jo 

e 


2(£o)^WpO,  +  eo  v,E\f 


(21) 


Since  space  charge  is  an  important  part  of  the  multipactor  phenomenon,  the  previ¬ 
ously  discussed  analysis  includes  these  effects.  According  to  Valfells,  et  al.  [16],  the 
main  authority  dealing  with  space  charge  effects,  the  power  deposited  can  vary  by 
up  to  a  factor  of  4  compared  to  when  space  charge  is  neglected.  Many  early  research 
efforts,  such  as  Kishek,  et  al.  [12]  ignored  space  charge  effects;  however,  all  recent 
research  includes  them  in  their  analysis,  since  the  effects  have  been  experimentally 
and  computationally  demonstrated.  The  effects  of  incorporating  space  charge  are  ev¬ 
ident  when  looking  at  Equation  16,  since  the  number  density,  n,  now  varies  spatially 
in  the  ^-direction.  This  means  that  Enc  also  varies  spatially,  hence  the  potential 
experienced  by  an  electron  is  no  longer  constant. 
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While  it  has  been  determined  that  space  charge  is  the  dominant  factor  in  multi¬ 
pactor  saturation,  the  term  ’’saturation”  may  be  misleading.  In  fact,  multipactor  is  a 
time-dependent  phenomenon,  as  shown  by  Kim  and  Verboncoeur  [9].  Since  all  work 
prior  to  this  discovery  relied  heavily  on  time-averaged  theoretical  developments,  the 
work  was  valid  for  that  assumption  and  can  be  used  as  a  guide  in  the  time-varying 
case,  as  will  be  shown  in  Section  2.2.4.  Kim  and  Verboncoeur  [9]  used  simulation 
methods  to  show  the  amplitude  of  the  electric  field  components  and  the  number  of 
electrons  vary  periodically  over  time.  This  has  also  been  observed  in  execution  of 
simulations  [7]  and  a  detailed  comparison  will  be  made  in  Chapter  III. 

2.2.4  Interpretation  and  Characterization  of  Multipactor. 

One  of  the  most  intuitive  and  useful  ways  to  characterize  multipactor  is  the  sus¬ 
ceptibility  curve.  Plots  such  as  the  one  given  in  Figure  6  are  commonly  referred  to  as 
susceptibility  curves.  The  lines  on  the  plot  represent  the  values  of  the  RF  field  (Erf) 
versus  the  restoring  field  (Erc)  f°r  which  the  SEY  is  one.  This  particular  suscep¬ 
tibility  plot  was  created  for  several  different  values  of  SmaX) 0,  which  would  represent 
different  types  of  materials.  The  axes  are  scaled  to  /  and  £max o>  so  that  the  suscep¬ 
tibility  curve  varies  only  for  values  of  the  secondary  emission  coefficient;  therefore, 
it  is  called  ’’universal.”  In  the  region  between  the  upper  and  lower  curves,  the  SEY 
is  greater  than  one;  therefore  multipactor  occurs,  resulting  in  an  electron  avalanche. 
When  the  corresponding  values  of  the  RF  field  and  the  restoring  field  fall  in  the  region 
above  the  upper  limit  or  below  the  lower  limit,  the  SEY  is  less  than  one  and  results 
in  total  loss  of  particles  back  to  the  dielectric  surface. 

It  is  worth  mentioning  the  conditions  under  which  these  susceptibility  curves  are 
generated  through  simulations.  For  all  cases  in  which  space  charge  is  ignored,  the 
value  for  Erc  is  assumed  to  remain  constant  and  multiple  calculations  are  made  with 
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Figure  6.  Reproduction  of  Kishek’s  universal  susceptibility  curves  [12]  for  the  simple 
(no  space  charge)  case.  The  values  are  normalized  to  the  frequency  and  Smaxo,  hence, 
the  term  ’’universal”. 

different  values  of  Erf-  The  pair  of  Er>c  /  Erf  values  is  then  characterized  based  on 
a  binary  system  (i.e. ,  either  multipactor  growth  occured  or  it  did  not).  Then,  Er>c 
is  changed  and  run  against  the  same  range  of  Erf  and  the  process  repeated  until 
the  multipactor  susceptibility  regions  are  well-defined.  This  is  a  valid  method,  under 
the  following  two  assumptions:  that  the  hoptime  of  an  electron  is  short  compared 
with  the  period  of  Erf  and  Er>c  does  not  change  spatially  (as  the  electron  moves 
away  from  the  surface).  The  latter  is  valid  based  solely  on  the  non-space  charge 
assumption. 

Analytical  expressions  have  been  developed  to  describe  these  types  of  susceptibil¬ 
ity  curves  by  Kishek  and  Lau  [11],  as  well  as  Semenov,  et  al.  [14],  While  the  two 
forms  of  the  boundaries  are  derived  in  different  ways  and  yield  seemingly  different 
forms,  there  is  good  agreement  on  the  boundaries.  The  two  formulas  for  multipactor 
boundaries  are  examined  below. 
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Kishek  and  Lau  [11]  give: 


eE 


RF 


LU 


y/m£, 


max  0 
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2  £c/£t 


max  0 


1  —  cos 


l  &u)y/m£or, 
\  eEflc 


(22) 


Where  £om  is  the  peak  value  for  the  emission  energy  distribution  of  the  electron, 
Erf  is  simply  the  value  for  the  RF  field  at  the  multipactor  boundary  (either  upper  or 
lower)  and  £c  is  the  crossover  energy  associated  with  the  same  boundary.  Sazontov 
and  Semenov  give 


2Wi/m  <  VRF  <  2 lb 2 /rn  for  weak  DC  fields 

( 2Wi/mVF )  V'yxj  <  VRF  <  (2W-2/jmVF)VFC  for  strong  DC  fields 


(23) 


Where  V  =  and  W\  and  W 2  are  the  crossover  energies.  In  order  to  compare  the 
formulas,  a  few  simplifying  assumptions  must  be  made.  For  example,  in  the  case  of 
a  strong  DC  field,  equation  22  becomes: 


eERF 

LO  y/  m£maxO 


2  £c/£, 


max  0 


1  f  4u\/m£om  A  2 

2  y  eEDC  ) 


(24) 


which  can  be  simplified  to  (after  noting  that  £om  =  kT)  : 


ErF  — 


(25) 


Similarly,  noting  that  Vt  =  y  ~~  for  strong  DC  fields,  one  side  of  the  inequality  in 
23  becomes  (using  Wc  to  denote  the  crossover  energy): 


ErF  — 


Edc  Iwc 
"\7§V  ~kT 


(26) 


Thus,  we  see  that  25  and  26  are  nearly  equal,  within  a  small  constant  factor.  Re- 
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gardless  of  this  small  constant  factor,  it  is  obvious  that  both  confirm  the  dependence 
of  Erf  on  Edc ,  £1,  and  kT. 

In  his  interpretation  of  the  susceptibility  curves,  Ang  [1]  notes  that,  in  contrast 
to  the  cavity  (waveguide)  case  for  which  beam  loading  was  the  primary  saturation 
mechanism,  the  primary  mechanism  leading  to  saturation  of  electrons  in  the  case  of 
a  single  dielectric  surface  is  the  growth  of  the  Edc  held  to  a  value  for  which  the  SEY 
is  unity.  He  also  notes  that  the  amount  of  energy  the  multipactoring  electrons  absorb 
is  a  small,  constant  percentage  of  the  total  RF  energy  (~  1%).  Although  this  may 
seem  insignificant,  in  the  case  of  HPM  systems  for  which  the  power  output  is  on  the 
order  of  gigawatts,  one  percent  of  the  RF  energy  would  indeed  be  enough  to  cause 
failure  in  the  dielectric  window. 

For  the  case  in  which  space  charge  effects  are  included,  the  susceptibility  curves 
(and  interpretation  thereof)  change  dramatically.  Kim  and  Verboncoeur  [9]  reveal 
that  a  more  accurate  representation  of  the  dynamics  of  a  multipactoring  system  is 
given  by  Figure  7  and  is  referred  to  as  the  ’’bowtie”  diagram  for  a  given  system. 

In  this  case,  the  dotted  lines  are  from  the  case  in  which  space  charge  was  ignored 
and  serve  as  ’’guidelines”  to  denote  the  region  in  which  multipactor  occurs.  These 
guidelines  have  been  reflected  around  the  Erf  =0  line  to  show  the  complete  variation 
of  the  RF  held,  instead  of  just  the  absolute  value.  The  plot  shows  the  tra  jectory  of  the 
system  in  susceptibility  space  ( Erf  versus  Edc  )  over  time.  The  t  —  0  occurs  at  the 
point  {Erc  =  0,  Erf  =  0)  and  as  time  progresses  in  the  simulation  the  trajectory 
of  the  system  follows  the  blue  curve.  The  regions  in  which  multipactor  occur  are 
between  the  top  two  red  dotted  lines  and  the  bottom  two  dotted  lines.  Therefore, 
the  values  of  ERF  and  EDC  never  reach  a  constant  steady  state,  as  was  the  case 
when  space  charge  effects  were  neglected.  Rather,  they  reach  a  periodic  steady  state, 
moving  in  and  out  of  the  multipactor  region.  As  would  be  expected  from  this  result, 
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the  particle  count  also  varies  periodically  (at  a  frequency  of  twice  the  period  of  ERF). 


Figure  7.  Reproduction  of  Kim  and  Verboncoeur’s  [9]  periodic  susceptibility  curves 
or  ’’bowtie”  diagrams.  Since  this  curve  is  not  normalized,  it  varies  for  each  set  of 
parameters  (ERf  ,  /,  8maxo,  £max o)-  The  parameters  for  this  plot  are:  Erfo  =  3  MV/m, 
/  =  1  GHz,  Smax o  =  2,  £maxo  =  400  eV.  The  dotted  lines  are  ’’guidlines”  from  the  case  in 
which  space  charge  effects  were  ignored. 


2.2.4. 1  Simulation. 

Since  experiments  are  extremely  difficult  and  expensive  to  conduct,  simulations 
are  the  most  practical  and  widely  used  method  of  examining  multipactor  in  HPM 
(and  other)  devices.  A  specific  type  of  simulation  code,  known  as  Particle-in-Cell 
(PIC)  code,  is  particularly  useful,  since  it  makes  very  few  physical  assumptions  and 
lends  itself  to  a  quantitative  analysis.  In  general,  PIC  codes  follow  the  mathematical 
flow  given  in  Figure  8,  obtained  from  Barker  [3]. 

The  first  major  simplification  PIC  code  introduces  is  to  represent  a  large  number 
of  particles  by  a  single  macropariticle,  which  maintains  the  same  charge  to  mass  ratio 
as  the  original  paritcle.  In  this  way,  the  number  of  particles  which  the  PIC  code 
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Figure  8.  The  general  mathematical  flow  of  PIC  codes,  as  given  by  Birdsall  [5]. 

must  track  and  account  for  is  inherently  reduced.  In  Figure  8,  the  index  i  is  used 
to  reference  quantities  associated  with  individual  particles  and  the  index  j  is  used 
to  reference  the  nodes  of  the  spatial  grid.  Therefore,  the  code  starts  by  solving  the 
equations  of  motion  for  each  particle.  Then,  since  the  position  of  each  particle  is 
known  exactly,  the  charge  of  each  particle  is  weighted  to  the  surrounding  corners  of 
the  grid.  Once  this  weighting  process  is  complete  for  every  particle,  the  total  charge 
contained  in  the  system  can  be  represented  by  a  charge  density  at  each  node,  pn  and 
a  current  density,  Jj,  between  the  nodes.  At  this  point,  the  Fast  Fourier  Transform 
(FFT)  is  used  to  solve  Maxwell’s  equations  (specifically,  Ampere’s  law  and  Faraday’s 
law)  for  the  corresponding  electric  and  magnetic  fields.  These  fields  can  then  be 
weighted  back  to  the  original  particle  locations  and  the  particles  moved  according 
to  the  force  exerted  on  them  by  the  electric  field.  Finally,  the  cycle  is  repeated 
as  the  simulation  continues.  The  Improved  Concurrent  Electromagnetic  Particle  in- 
Ccll  (ICEPIC)  [4]  code  is  remarkable  in  it’s  ability  to  self-consistently  calculate  and 
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model  all  the  electromagnetic  effects  within  a  system,  both  internally  generated  by 
the  particle  interactions  and  externally  imposed.  In  the  course  of  this  research,  two 
dimensional  simulations  (with  three  dimensional  velocity  vectors)  are  run  in  order  to 
maximize  computational  efficiency. 

There  are  a  few  restrictions  and  conditions  that  apply  to  general  PIC  code.  The 
Courant  condition  governs  the  time  step  (dt)  associated  with  the  simulation.  Accord¬ 
ing  to  this  condition,  the  time  step  depends  on  the  size  of  the  cell  and  the  relevant 
velocity.  In  the  case  in  which  space  charge  is  ignored,  this  velocity  is  Vt,  the  thermal 
velocity.  However,  in  the  space  charge  case,  the  relevant  velocity  is  the  speed  of  light, 
Co,  which  is  several  orders  of  magnitude  higher.  This  means  that  the  grid  must  be 
considerably  more  refined  (i.e.  smaller  and  a  greater  number  of  cells)  for  the  space 
charge  case.  Additionally,  it  is  important,  in  order  to  take  advantage  of  ICEPIC’s 
concurrent  electromagnetic  capabilities,  the  cell  dimension  should  be  on  the  order  of 
the  Debye  length.  Multipactor  is  also  a  statistical  phenomenon,  therefore,  for  simu¬ 
lations  to  be  statistically  valid,  each  cell  in  the  interaction  space  should  contain  5-10 
macroparticles. 

Computational  PIC  tools,  such  as  ICEPIC  and  XPDP  [19],  in  conjunction  with 
Monte-Carlo  simulations  are  the  predominant  modeling  codes  used  to  verify  analytical 
developments.  While  the  development  of  the  analytical  theory  is  well  documented  by 
Kishek,  et  al.  [10],  Ang,  et  al.  [1],  Valfclls,  et  ah  [16],  and  Kim  and  Verboncoeur  [9], 
the  simulation  parameters  were  not  as  thoroughly  documented  throughout  the  same 
literature.  However,  in  general,  the  results  agree  well  with  the  theory.  According  to 
Ang,  et  ah  [1]  ,  typical  parameters  for  simulations  are: 

dmax o  =  3  ,Smax o  =  420eV,  Som  =  2.1eV,  Erfq  =  3MV/m 
/  =  2.45GHz,  Edco  =  lOV/m,  Z0  =  377H 
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Using  a  Monte-Carlo  simulation,  Ang,  et  al.  [1]  were  able  to  demonstrate  that,  after 
an  initial  (very  short)  transient  period,  parameters  such  as  n  (t),  E^c,  5  (t),  Si,  t  and 
Pdep/ Prf  settled  into  periodic  steady  state  values  near  those  that  were  predicted  by 
analytical  theory.  This  is  consistent  with  later  research  demonstrated  by  Kim  and 
Verboncoeur  [9]  and  Fichtl  [7]. 

Valfclls,  et  al.  [16]  were  also  able  to  show  good  correlation  with  theoretical  results, 
using  the  XPDP1  code.  Fichtl  [7]  updated  ICEPIC  to  include  secondary  emission 
effects,  verifying  the  code  by  comparing  it  to  theory.  Discoveries  made  during  the 
course  of  this  research  will  be  applied  to  update  and  correct  Fichtl’s  ICEPIC  SEE 
code. 

2.2.5  Mitigation  Methods. 

Mitigation  of  multipactor  in  HPM  systems  is  one  of  great  current  interest  among  the 
leading  experts  in  the  field.  Investigation  of  two  promising  mitigation  methods  is 
the  primary  aim  of  this  paper:  depositing  a  thin  coating  of  low-SEY  material  on  the 
dielectric  window  and  geometrical  modification  of  the  dielectric  window. 

2. 2. 5.1  Thin  Coatings  of  Low-SEY  Materials. 

Since  the  materials  used  to  construct  the  dielectric  window  in  a  HPM  system 
typically  have  Smax o  ~  2  —  3,  one  possible  method  of  multipactor  mitigation  is  to 
reduce  the  SEY  by  deposition  of  a  thin,  durable  and  low-SEY  material  on  the  surface 
of  the  dielectric.  This  accomplishes  two  objectives:  the  previously  mentioned  lowering 
°f  5 max o  (arid,  consequently,  Smax )  and  the  ”  squeezing  in1'  of  the  crossover  energies, 
so  that  the  range  of  impact  energies  for  which  the  SEY  is  greater  than  1  is  much 
smaller.  Both  of  these  effects  are  evident  from  Figure  9.  Desirable  properties  for  such 
a  material  include:  1)  low  SEY,  2)  good  thermal  stability  and  3)  minimum  attenuation 
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of  Erf  ,  so  that  the  him  does  not  become  yet  another  significant  loss  mechanism. 
According  to  Castanada,  et  ah  [6]  and  Kennedy,  et  al.  [8],  the  SEY  of  Titanium 
Nitride  (Ti:N)  is  much  lower  than  2  and  can  even  be  deposited  and  processed  in  a 
manner  that  brings  it  near  1.  Of  course,  this  is  an  ideal  situation,  but  such  processing 
is  very  difficult  and  precise  work;  therefore,  it  must  be  accomplished  carefully  in 
pristine,  vacuum  conditions.  If  the  real  secondary  emission  curve  for  such  a  material 
were  known,  it  would  be  easy  to  simulate.  However,  secondary  emission  curves  are 
also  very  difficult  to  measure  and  require  an  elaborate  apparatus,  so  we  can  currently 
only  estimate  the  parameters  of  materials.  As  can  be  seen  from  Equations  1-5,  the 
only  parameters  necessary  to  generate  a  secondary  emission  curve  for  a  given  material 
are  5maX:0  and  £maxo  (since  (  and  ks  are  usually  assumed  to  remain  constant).  Then, 
using  ICEPIC,  simulations  can  be  run  to  determine  the  susceptibility  to  multipactor 
of  the  new,  coated  system. 

2. 2. 5. 2  Geometrical  Modification  of  the  Waveguide/ Window  Inter¬ 
face. 

Modifying  the  geometry  of  the  interface  of  the  dielectric  window  with  the  waveg¬ 
uide  means,  in  the  simplest  case,  introducing  a  non-zero  angle  of  incidence  between 
the  window  and  Erf  •  This  is  equivalent  to  the  window  no  longer  being  perpendic¬ 
ular  to  the  waveguide,  but  offset  by  the  angle  i/j.  With  regard  to  RF  field  oblique 
orientation,  Valfclls,  et  al.  [15]  begin  with  the  fact  that  the  electron  flight  time  is  the 
most  important  factor  in  determining  the  susceptibility  of  a  system  to  multipactor 
damage.  In  general,  the  longer  an  electron  is  in  flight,  the  more  energy  it  gains  from 
the  RF  field  and  the  more  energy  it  possesses  when  it  returns  to  the  surface.  An 
oblique  RF  field  angle  of  incidence,  introduces  a  component  of  Erf  in  the  direction 
of  Edc  and  can  either  enhance  or  degrade  the  time  of  flight  of  electrons,  thereby  in- 
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Secondary  Emission  Yield 


Figure  9.  A  comparision  of  secondary  emission  yield  curves.  The  only  parameter 
varied  in  this  plot  is  Smax o,  showing  the  drastic  difference  a  reduction  in  this  material 
parameter  can  make. 


creasing  or  decreasing  the  susceptibility  of  the  system  to  multipactor.  Such  a  system 
is  one  in  which  0  in  Figure  10  is  greater  than  zero.  According  to  Valfells,  et  al.  [15], 
for  0  <  ~  5°,  the  susceptibility  curve  is  widened;  however,  larger  values  of  0  result 
in  a  considerably  narrower  susceptibility  curve. 

Formulas  for  the  trajectories  and  flight  times  of  electrons  in  a  system  with  an 
obliquely  incident  RF  field  are  given  by  Valfells,  et  al.  [15]: 

vx  =  v0  sin (0)  -  EDCt  +  Erfq  sin(0)  (cos (t  +  9)  -  cos(0)) 
vy  =  Vo  cos(0)  +  Erfo  cos(0)  (cos(£  +  6)  —  cos{6 )) 

0  =  (v0  sin(0)  -  Erfo  sin(0)  cos(6>))  T  -  T 2 

+ERFos'm(^)  (sin(T  +  6)  -  sin(61)) 

Formulas  27  -  29  can  be  used  to  calculate  the  energy  of  an  electron  returning  to 
the  surface  using  the  elementary  relationship  1/2  (n/  +  v^j  [Note:  The  mass  of  the 


(27) 

(28) 
(29) 
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Figure  10.  The  geometry  for  single  surface  multipactor.  This  is  the  same  as  Figure 
2,  except  an  angle  of  incidence,  %/>,  has  now  been  introduced  between  the  surface  and 
Erf  •  This  is  equivalent  to  placing  the  window  at  some  angle  if>  off  of  normal  to  the 
waveguide. 


electron  has  been  normalized  out  of  this  set  of  equations.]  Valfclls,  et  ah  [15],  use 
simulations  to  obtain  the  susceptibilty  curves  for  the  oblique  incidence  case.  These 
simulations  revealed  that  for  values  of  ip  <  5  —  10°,  the  susceptibility  curve  is  broad¬ 
ened,  but  for  larger  values  of  i/j,  the  region  of  susceptibility  to  multipactor  is  drastically 
reduced.  The  reason  for  the  broadening  at  small  angles  of  obliqueness  is  related  to 
the  slope  of  the  secondary  emission  curve  at  the  lower  boundary.  Figure  11  shows 
the  SEY  plotted  as  a  function  of  the  restoring  held  (. Edc  ),  rather  than  the  impact 
energy.  The  higher  restoring  held  is  related  to  the  lower  boundary  on  the  suscepti¬ 
bility  curve  (which  is  the  hrst  crossover  energy,  E\.  When  a  component  of  Erf  exists 
in  the  direction  of  the  Edc  ,  the  restoring  held  can  now  be  either  slightly  smaller  or 
slightly  larger,  depending  on  the  phase  of  the  RF  held.  Because  the  slope  is  decreas- 
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ing  near  this  point,  denoted  by  the  lower  boundary  label  in  Figure  11,  the  average 
secondary  emission  due  to  this  small  change  in  the  restoring  field  would  be  higher 
than  for  the  restoring  field  by  itself.  For  larger  angles  of  obliqueness,  the  RF  field 
component  in  the  direction  of  the  restoring  field  generally  decreases  the  restoring  field 
by  a  greater  amount,  resulting  in  a  longer  flight  time  and  electrons  returning  to  the 
dielectric  surface  with  more  energy,  such  that  the  SEY  becomes  less  than  one.  This 
effect  dominates  over  the  previously  discussed  change  in  the  restoring  field.  Thus, 
the  susepetibility  region  will  be  significantly  diminished  for  larger  angles  of  oblique¬ 
ness.  An  analytical  method  which  utilizes  MATLAB  to  examine  the  change  in  the 
susceptibility  curves  is  described  in  detail  in  Chapter  Ill  and  the  results  compared 
with  Valfells’  in  Chapter  IV. 


Figure  11.  Secondary  emission  curve  as  given  by  Valfells,  et  al.  [15] 


III.  Computational  Models  and  Comparison  with  Theory 


This  chapter  describes  the  rationale  behind  the  selection  of  ICEPIC  input  pa¬ 
rameters.  including  assumptions  made  in  running  simulations.  This  is  followed  by  a 
detailed  description  of  basic  simulations  that  were  run  in  order  to  verify  our  system 
by  comparing  the  results  with  the  theory  described  Kishek,  et  al.  [10],  Ang,  et  al. 
[1],  Valfells,  et  al.  [16]  and  Kim  and  Verboncoeur  [16].  Also,  the  numerical  and 
simulation  models  for  potential  mitigation  methods  are  thoroughly  described. 

3.1  ICEPIC  Input  Parameters 

This  section  includes  a  detailed  discussion  and  explanation  of  the  input  parame¬ 
ters  used  in  simulating  the  cases  follows,  broken  up  into  functional  sections.  While  a 
nominal  understanding  of  the  operation  of  ICEPIC  is  necessary  to  completely  grasp 
the  complexity  of  the  system,  a  general  understanding  of  PIC  code  will  suffice  when 
reviewing  the  following  discussion.  ICEPIC  simulations  are  controlled  by  input  hies, 
which  contain  parameters  that  allow  a  wide  variety  of  physical  phenomena  to  be  sim¬ 
ulated.  This  section  will  describe  each  section  used  in  generating  the  simulations  used 
throughout  the  course  of  our  research,  explaining  the  purpose  of  the  most  important 
parameters.  As  a  matter  of  convention,  each  functional  section  in  the  input  hie  is 
denoted  by  a  title  enclosed  in  square  brackets  -  i.e,  the  hrst  section  is  the  [Defaults] 
section.  Also,  comments  are  denoted  by  a  semicolon,  such  that  everything  from  the 
semicolon  to  the  end  of  the  current  line  is  ignored  by  the  code.  When  running  an 
input  hie,  ICEPIC  considers  the  entire  hie  and  is  not,  for  the  most  part,  sensitive 
to  the  order  in  which  sections  are  entered.  The  complete  input  hies  for  both  the 
NSC  case  and  the  SC  case  can  be  found  in  Appendices  A  and  B.  For  a  complete 
explanation  of  every  available  parameter  and  section,  please  see  the  ICEPIC  User’s 
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Manual  [4]. 


3.1.1  Defaults  and  Variables. 

The  [Defaults]  and  [Variables]  sections  allow  a  user  the  ability  to  define 
variables  for  easy  reference  and  usage  throughout  the  rest  of  the  input  file.  The 
difference  between  the  two  sections  is  that  variables  placed  in  the  [Defaults]  section 
are  eligible  for  parametric  surveys.  In  this  case,  a  separate  script  can  be  used  to 
generate  many  input  files  containing  a  survey  of  values  for  a  particular  variable.  As 
an  example,  the  input  files  given  in  Appendix  A  contains  the  ENERGY  variable,  which 
is  later  used  to  assign  the  value  of  Erf  (for  the  NSC  case).  In  building  the  surveys 
that  were  used  to  create  the  lower  multipactor  boundary,  a  script  was  run  to  generate 
many  input  files  containing  a  range  of  values  of  for  both  Erf  and  E^c  ■ 

3.1.2  [Time]  Section. 

The  [Time]  section  contains  parameters  that  pertain  to  the  simulation  time,  in¬ 
cluding: 

•  dt  -  This  is  the  fundamental  time  step,  of  the  simulation. 

•  stepunax  -  This  represents  the  total  time  for  the  simulation,  in  units  of  number 
of  dt  steps. 

•  c  our  ant  .value  -  This  is  a  constant  value  which  determines  what  the  time  step 
should  be  set  to  should  the  Courant  condition  not  be  satisfied. 

•  kill.if  _below_parts  -  Once  the  number  of  particles  in  the  system  reaches  this 
minimum  value,  the  simulation  is  terminated. 

•  kill.if  _above_parts  -  Once  the  number  of  particles  in  the  system  reaches  this 
maximum  value,  the  simulation  is  terminated. 
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3.1.3  System  Geometry. 


The  system  geometry  is  defined  in  many  different  sections  of  the  input  file,  specif¬ 
ically: 

•  [Cartesian]  -  Setting  ig_dir  to  (0,1,0)  causes  the  code  to  physically  ignore 
the  y  direction,  making  this  a  two-dimensional  survey.  In  the  NSC  case,  E^c  is 
set  by  the  default _e  parameter. 

•  [XGridl] ,  [YGridl] ,  [ZGridl]  -  These  three  parameters  define  the  computa¬ 
tional  grid  for  the  simulation,  including  the  physical  extent,  the  size  of  the  cells, 
dx,  dy  and  dz,  and  the  number  of  cells  in  each  direction,  Nx,  Ny  and  Nz. 

•  [Symmetry]  -  Sets  the  periodic  boundary  conditions  for  the  system. 

•  [ShapeN]  -  These  sections  define  physical  shapes  in  the  computational  domain. 
The  initial  domain  is  considered  to  be  a  metallic  cube  and  the  physical  shapes 
are  ”  carved”  out  of  this  cube  using  a  number  of  [ShapeN]  sections.  The  sec¬ 
tions  are  numbered  during  the  pre-simulation  processing  script.  The  physical 
properties  of  the  shape  are  determined  by  the  material  parameter,  which  may 
be  more  specifically  defined  in  the  [MaterialN]  section. 

•  [MaterialN]  -  The  physical  properties  of  the  shapes  may  be  defined  and  refer¬ 
enced  by  assigning  a  name  to  any  number  of  materials.  Similar  to  the  [ShapeN] 
sections,  these  sections  are  numbered  during  the  pre-simulation  processing. 

•  [BoundN]  -  Boundary  conditions  may  be  set  using  these  sections.  For  the  pur¬ 
poses  of  this  research,  a  Perfectly  Matched  Layer  (PML)  boundary  condition  is 
employed. 

Figure  12  shows  the  geometrical  configuration  of  the  system  used  during  these  simula¬ 
tions,  which  is  generated  by  the  above  sections.  This  is  very  similar  to  the  geometries 
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used  in  previous  work  by  Fichtl  [7]  and  Verboncoeur,  et  al.  [9].  In  this  configuration, 
the  pyhsical  size  of  the  system  scales  with  the  wavelength  in  any  direction.  The  di¬ 
electric  material  is  given  a  default  eo  of  1.001,  so  that,  in  the  electromagnetic  sense, 
it  simulates  a  nearly  opaque  material,  allowing  the  RF  energy  to  propagate  with  no 
significant  attenuation.  This  is,  naturally,  a  realistic  assumption  for  HPM  systems. 
The  purpose  of  the  dielectric  on  the  left  side  of  the  system  is  to  shield  the  plane  wave 
emitter  from  any  particles  that  may  traverse  the  length  of  the  system.  For  the  pa¬ 
rameters  given  in  Section  2. 2. 4.1,  typical  values  for  the  SC  case  would  be  (for  Nx  =  3 
and  Nz  —  500): 
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In  order  to  fulfill  the  condition  that  there  must  be  many  particles  contained  in  a 
Debye  sphere  when  considering  individual  particle  interactions,  dz  for  the  SC  case  is 
four  times  smaller  than  dz  for  the  NSC  case,  therefore  the  dimensions  for  the  NSC 
case  are  four  times  larger  than  those  in  the  SC  case. 

A  PML  enforces  the  electromagnetic  boundary  conditions  behind  the  dielectric 
and  prevents  any  backscattered  waves  from  affecting  the  simulation  by  introducing 
layers  of  successively  higher  conductivity.  The  system  is  symmetrical  (or  periodic)  in 
the  ^-direction,  so  that  any  particle  leaving  the  system  in  on  one  side  of  the  system 
at  xrnax  reenters  the  system  at  xmin  and  vice-versa. 
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Ax  = 
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Figure  12.  The  geometry  used  for  ICEPIC  simulations.  Nx  and  Nz  are  the  number  of 
cells  in  the  respective  direction,  dz  =  A/3000  is  the  width  of  a  cell.  The  particle  emitter, 
planewave  emitter  and  planewave  absorber  take  up  no  physical  space. 

3.1.4  Plane  Wave  Emission. 

Erf  is  simulated  by  a  plane  wave  that  extends  across  the  system  in  the  ^-direction, 
originates  from  a  plane  behind  the  left  hand  dielectric  and  propagates  in  the  positive  ,0- 
direction,  as  shown  in  Figure  12.  In  the  input  hie,  this  is  created  by  the  [PlanewaveN] 
section.  The  ICEPIC  code  actually  creates  a  plane  wave  emitter  on  one  side  of  the 
system  and  a  plane  wave  absorber  on  the  other  side  of  the  system,  for  computational 
stability.  This  is  reflected  in  the  definition  of  the  shape  parameter  in  this  section. 

The  magnitude  of  Erf  is  determined  by  the  EO  parameter,  which  is  set  by  the 
energy  variable,  which  is  assigned  a  default  value  in  the  [Defaults]  section  by 
ENERGY.  The  frequency  of  the  system  is  taken  to  be  that  of  the  RF  held.  For  the  NSC 
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case,  the  restoring  field  (E^c )  is  set  in  the  [Cartesian]  section  by  the  ez  variable, 
which  is  also  assigned  a  default  value  in  the  [Defaults]  section  by  EDC.  This  same 
methodology  is  used  to  create  a  time-varying  E-field  in  the  case  of  the  geometrical 
modification  mitigation  method  (the  specific  implementation  is  discussed  in  Section 
3.3). 

3.1.5  Particle  Emission. 

Particle  emission  is  generated  in  the  ICEPIC  input  file  in  the  [ParticlesN]  sec¬ 
tions.  Two  main  variants  of  this  section  are  used  through  the  course  of  this  research: 

•  Primary  particle  emitter  -  This  section  provides  seed  particles  to  the  system 
and  is  defined  by  setting  method  to  BEAM.  The  parameters  that  control  the 
beam  emission  are: 

—  temp  -  Controls  the  initial  temperature  of  the  injected  particles. 

—  inject.interval  -  Determines  the  injection  frequency  of  the  particles  and 
is  given  in  units  of  dt. 

—  current  -  One  of  the  main  factors  in  determining  the  macroparticle  weight¬ 
ing  -  which  is,  in  this  case,  is  fully  described  by:  where  inj  ect 

is  the  number  of  particles  emitted  per  emitting  face  (where  each  Nx  is  an 
emitting  face)  per  time  step  and,  for  this  case,  is  1. 

—  tstop  -  Determines  the  time  at  which  particle  injection  ceases,  in  this  case 
particles  are  injected  for  half  an  RF  period. 

—  The  other  parameters  contained  in  this  section  are  used  to  smooth  out  the 
emission  waveform  and  probably  may  be  omitted. 

•  Secondary  particle  emitter  -  This  section,  specified  by  setting  method  to  SECONDARY, 
defines  the  secondary  emission  properties  of  the  dielectric  face,  per  the  sec- 


35 


ondary  emission  equations  given  in  Section  2.2.2.  The  parameters  that  control 
the  secondary  emission  properties  are: 

—  sec_coef  f  icient  -  Corresponds  to  5max o. 

—  sec_energy  -  Represents  the  average  emission  energy  of  secondary  parti¬ 
cles. 

—  sec_ref  1  -  The  fraction  of  secondary  electrons  that  are  reflected  from  the 
surface  of  the  dielectric. 

—  sec_scat  -  The  fraction  of  secondary  electrons  that  are  scattered  from  the 
surface  of  the  dielectric. 

—  sec_energynnax_yield  -  Corresponds  to  £maxo- 

—  sec_ks  -  The  smoothing  constant,  related  to  the  surface  smoothness  of  the 
material. 

In  each  section,  the  properties  of  the  particles  are  defined  by  the  mass  and  q 
(charge)  parameters.  The  dir  parameter  defines  the  direction  in  which  the  particles 
are  emitted.  Also,  the  shape  parameter  defines  the  plane  face  from  which  the  particles 
will  be  emitted. 

The  particle  emitter  and  secondary  emitter  are  created  by  very  narrow  (about  half 
a  cell  wide),  non-physical  faces  on  the  surface  of  the  right  hand  dielectric.  The  primary 
particle  emitter  section  is  used  to  provide  a  small  amount  of  seed  current  to  the 
system.  [NOTE:  in  ICEPIC,  the  current  parameter  is  one  of  the  parameters  used  to 
determine  the  macroparticle  weighting,  that  is,  how  much  charge  each  macroparticle 
carries.]  This  current  is  large  enough  to  provide  seed  electrons  to  the  system,  but  not 
so  large  as  to  affect  whether  or  not  the  system  is  operating  in  a  multipactoring  region 
or  to  produce  an  electron  density  that  is  too  high  (such  that  a  large  amount  of  power 
is  reflected  back  towards  the  source).  Also,  the  current  is  only  injected  for  half  an 
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RF  period,  so  that  it  will  not  interfere  with  the  amount  of  energy  measured  in  the 
particles  or  the  power  of  the  system.  The  particles  are  injected  with  a  Maxwellian 
distribution  of  velocities,  given  by: 


/  M  = 


2  n0 

V^vt 


exp 


where  vt  =  \J  therefore  the  energy  of  the  seed  electrons  is  specified  by  a  particular 
temperature,  T,  which  is  the  parameter  temp  in  the  primary  particle  emitter  section 
of  the  ICEPIC  input  hie.  The  length  of  time  for  the  seed  period  is  one-half  of  an 
RF  period,  so  as  to  inject  particles  at  many  phases  of  the  RF  held.  In  this  case,  the 
particles  are  injected  with  a  frequency  of  30  times  per  RF  period. 

Several  simulations  were  run  to  determine  the  most  effecient  way  to  control  particle 
weighting  and,  effectively,  the  run  time  of  the  simulation.  It  was  found  that,  for  a 
constant  emitting  area,  that  maintaining  the  default  injection  rate  of  1  macroparticle 
per  nx  per  dt  and  increasing  the  injection  current  results  in  a  proportional  increase  in 
the  macroparticle  charge.  This  further  results  in  a  decreased  number  of  particles  in  the 
system,  higher  computational  efficiency  and  strikingly  comparable  results.  This  does 
not  cause  any  distortion  of  the  results  because  the  macroparticle  charge  to  electron 
charge  ratio  is  always  accounted  for  in  the  post-simulation  analysis  of  data.  However, 
for  the  case  in  which  Erf  =  5.5  MV/m,  this  resulted  in  a  run-time  savings  of  almost 
40  hours,  since  ICEPIC  run  time  increases  nearly  exponentially  with  the  amount  of 
particles  in  the  system.  The  most  important  factor  to  consider  is  to  make  sure  that 
the  critical  electron  density  is  not  exceeded,  such  that  the  reflected  power  is  not  too 
high.  The  critical  electron  density  is  the  high  electron  density  at  which  total  reflection 
of  incoming  power  occurs  and  is  discussed  in  further  detail  in  Section  4.1.2.  Also, 
enough  macroparticles  must  be  present,  so  as  to  produce  a  statistically  valid  answer. 
Typically,  this  is  on  the  order  of  a  few  macroparticles  per  cell  over  the  interaction 
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space.  Therefore,  determining  the  appropriate  macroparticle  weighting  is  a  delicate 
balancing  act  and  must  be  carefully  considered  when  running  PIC  simulations. 

The  secondary  emitter  surface  models  the  secondary  emission,  as  described  by 
Vaughn  [17] .  The  energy  of  the  emitted  secondaries  is  also  determined  by  a  Maxwellian 
distribution  and  is  specified  by  the  parameters  in  the  [ParticlesN]  section  denoted 
by  method=SECONDARY.  Scattering  and  reflection  effects  are  also  included  for  all  sim¬ 
ulation  and  analyses,  per  Fichtl  [7],  except  where  otherwise  specifically  noted. 

3.1.6  Dump  Parameters. 

In  ICEPIC,  many  different  quantities  may  be  written  to  output  hies,  or  ’’dumped” 
throughout  the  course  of  a  given  simulation,  using  the  [DumpN]  sections.  The  quan¬ 
tities  may  be  dumped  over  the  whole  simulation  or  during  any  specified  part  of  it. 
Specific  dumps  are  created  by  setting  method  to  the  dump  name.  As  with  other 
sections,  the  [DumpN]  sections  are  numbered  during  the  pre-simulation  processing 
scripts.  While  some  of  the  dump  sections  in  the  sample  input  hies  are  commented 
out,  the  most  useful  sections  are  described  below: 

•  FIELD  -  This  dumps  the  vector  electric  and  magnetic  held  quantities  at  every 
node  in  the  specihed  volume.  In  the  course  of  this  research,  the  RF  held  taken  to 
be  the  x-directed  component  and  the  restoring  held  is  taken  to  be  the  ^-directed 
component. 

•  POWER  -  The  Poynting  hux  through  the  plane  specihed  by  the  dumpwalue  pa¬ 
rameter  is  dumped.  This  represents  a  power  density  and  can  easily  be  converted 
to  power  using  the  system  geometry. 

•  PART  -  This  dumps  out  the  position,  velocity  and  charge  for  every  particle  in  the 
specihed  volume.  Since  the  number  of  particles  in  a  simulation  can  become  very 
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large,  these  dump  files  can  become  quite  large  and  should  be  used  sparingly. 
The  PART  dump  also  allows  for  only  information  about  particles  collected  at 
any  boundary  in  the  system  to  be  written,  which  drastically  reduces  the  size 
of  the  files  and  allows  calculation  of  the  power  deposited  on  the  surface  of  the 
dielectric  window,  as  desscribed  in  Section  3.2.3. 

Every  dump  parameter  may  be  controlled  by  a  common  set  of  parameters,  including: 

•  dump.int  -  Determines  how  frequently  dump  files  are  created. 

•  nstart  -  Sets  the  starting  time  for  the  dump  files,  in  units  of  (It. 

•  nstop  -  Sets  the  stop  time  for  the  dump  files,  in  units  of  dt. 

3.1.7  The  [Expert]  Section. 

The  [Expert]  section  contains  two  vital  parameters  used  through  in  performing 
these  simulations: 

•  requested_max_num_particles  -  Once  the  simulation  reaches  the  number  of 
particles  specified  by  this  parameter,  the  weighting  of  particles  emitted  by  the 
secondary  emission  code  will  be  increased  instead  of  emitting  more  than  one 
particle.  This  is  one  way  to  maximize  computational  efficiency,  but  makes  post¬ 
simulation  analysis  of  particles  much  more  difficult,  due  to  the  many  different 
weightings  that  may  be  present  in  the  system. 

•  self  .fields  -  This  parameter  determines  whether  or  not  space  charge  effects 
are  included.  If  set  to  0,  ICEPIC  simulates  the  NSC  case,  such  that  the  particles 
only  move  under  the  influence  of  external  fields.  If  set  to  1,  the  current  density 
of  the  particles  is  incorporated  in  the  solution  of  Maxwell’s  equations. 
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3.1.8  ICEPIC  Code  Updates. 


In  the  process  of  verifying  previous  results,  updates  to  ICEPIC  code  were  neces¬ 
sitated.  While  this  is  by  no  means  a  complete  listing  of  updates  made  to  the  code 
during  the  period  in  which  this  research  was  performed,  this  section  contains  a  short 
description  of  some  of  main  code  updates  that  were  required  through  the  course  of 
this  research. 

1.  DENSITY  Dump  -  It  was  discovered  that  this  particular  dump  contains  an  error 
and  produces  no  output  when  running  on  a  multi-processor  machine.  This  error 
did  not  affect  the  scope  of  our  research;  however,  the  DENSITY  dump  is  an  ex¬ 
tremely  useful  dump,  containing  (for  each  cell  and  dump  interval):  real  particle 
number  density,  velocity  information,  energy  information  and  the  number  of 
macroparticles. 

2.  POWER  dump  unit  conversion  -  A  factor  of  100  reduction  is  present  in  the  POWER 
dump  hie  values.  This  unit  conversion  comes  from  the  internal  dimension  for 
the  ignorable  direction  (in  this  case,  y )  being  set  to  1  cm  rather  than  1  m  in  a 
two  dimensional  simulation. 

3.  PART  dump  unit  conversion  -  A  factor  of  100  was  found  in  the  charge  column  of 
the  dump  hie,  when  collectecLpart=l.  The  values  reported  in  this  research 
have  been  compensated  for  this  factor. 

4.  Surface  Normal  Error  -  The  general  particle  emission  code  was  passing  an  in¬ 
correct  cell  normal  direction  to  the  secondary  electron  emission  code,  which  was 
then  using  that  value  to  determine  emission  and  reflection  angles  for  reflection 
and  scattering  of  secondary  electrons.  This  bug  was  hxed  in  the  version  of 
ICEPIC  used  for  this  work. 
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The  input  files  used  in  the  course  of  this  research  were  based  on  those  presented 
by  Fichtl  in  [7],  however,  several  parameters  in  Fichtl’s  input  files  were  no  longer 
supported  by  the  ICEPIC  code.  These  parameters  were: 

1.  symmetry,  symunode  in  the  [Symmetry]  section  -  These  legacy  parameters  were 
replaced  with  the  updated  periodic  parameters. 

2.  cell_mult  in  the  [Expert]  section  -  No  documentation  was  found  for  this 
parameter,  so  it  was  removed  from  the  code.  Apparently,  it  was  a  correction 
factor  for  a  known  issue  in  ICEPIC,  which  has  been  fixed. 

3.  [ShapeN]  sections  -  An  adjustable  buffer,  called  zbuff,  was  introduced  in  the 
SHAPE  sections  in  the  z-direction  on  both  sides  of  the  Plane  Wave  Emitter/Ab¬ 
sorber  (on  both  ends  of  the  system),  in  order  to  investigate  the  effects  of  the 
size  of  such  a  buffer  on  POWER  and  FIELD  dumps.  It  was  thought  that  several 
cells  on  either  side  of  the  Plane  Wave  Emitter/Absorber  allowed  more  room 
for  the  FIELD  and  POWER  dumps  and  left  less  potential  for  instabilities  caused 
by  the  proximity  of  the  boundaries  to  the  dump  planes.  However,  very  little 
difference  was  found  when  introducing  a  6 dz  wide  buffer  versus  the  1  dz  buffer 
used  by  Fichtl. 

3.2  Comparison  of  Results 

3.2.1  Non-Space  Charge  (NSC)  Case. 

The  NSC  case,  in  which  space  charge  effects  are  ignored,  is  the  simplest  case  and 
the  most  rational  starting  point  for  a  verification  of  multipactor  theory  in  ICEPIC. 
Most  of  the  theory  from  this  section  will  be  compared  with  Fichtl’s  [7]  results.  In 
this  case,  the  electromagnetic  contributions  of  particles  are  ignored  (meaning  that  the 
current  density  is  not  updated  in  Maxwell’s  equations)  and  only  the  fields  imposed  by 
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external  sources  are  considered.  Thus,  it  is  necessary  to  specify  a  constant  value  for 
Edc  in  order  to  facilitate  multipactor  growth.  Using  the  same  methodology  Fichtl  [7] 
uses  (which  is  described  in  Section  2.2.4),  the  lower  multipactor  boundaries  in  Figure 
13  were  reproduced.  The  upper  boundary  is  not  considered,  because,  according  to 
Barker  [3],  the  lower  boundary  is  associated  with  saturation.  In  this  case,  the  only 
metric  used  to  determine  multipactor  growth  was  whether  or  not  the  number  of 
particles  in  the  system  grows  significantly  beyond  the  initial  seed  particles.  In  order 
to  maximize  computational  efficiency,  the  simulations  are  terminated  if  they  reach  a 
maximum  number,  specified  by  the  kill_if  .above  parameter  in  the  [Time]  section 
of  the  input  file  in  .  This  number  is  well  above  the  number  of  particles  necessary 
to  qualify  a  particular  simulation  for  multipactor  growth.  One  notable  variation 
between  the  results  is  the  fact  that  the  boundaries  do  not  match  up  exactly.  This 
is  most  likely  due  to  updates  to  the  secondary  emission  code  within  ICEPIC  since 
Fichtl  performed  his  research.  However,  the  correlation  of  the  curves  is  satisfactory 
and  confirms  the  simulations  were  appropriately  configured.  This  process  could  be 
repeated  for  different  secondary  emission  coefficients,  however,  these  types  of  surveys 
are  time  consuming  and  it  is  assumed  that  this  result  implies  that  other  surveys  would 
demonstrate  similar  behavior. 

In  order  to  demonstrate  the  rapidly  evolving  nature  of  multipactor,  Figure  14 
shows  the  particle  growth  in  regions  near  the  lower  susceptibility  boundary.  From 
these  figures,  it  is  evident  that  multipactor  growth  occurs  very  quickly  -  usually 
within  a  few  nanoseconds  (where  the  period  of  the  RF  cycle  is  0.4  ns  and  a  hoptime 
may  be  on  the  order  of  1  ps).  These  points  at  which  these  plots  are  generated  are 
different  from  those  in  Figure  8  of  [7],  because  they  have  been  adjusted  to  correspond 
to  the  susceptibility  curves  generated  by  more  current  simulations  (i.e.  those  shown 
in  Figure  13). 
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Multipactor  Region  Boundaries 


Edc  [MV/m]  x(f/l  GHz)-1x(BraaI/400  eV)'1^ 

Figure  13.  Comparison  with  Fichtl’s  multipactor  boundary  curves.  The  labeled  points 
correspond  with  the  points  selected  for  Figure  14  and  are  in  reference  to  the  Rogers 
case  with  no  scattering  or  reflection  effects  (solid  green  line). 
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Edc  =  1-000  MV/m,  Erf  =  3.7  MV/m 


Edc  =  1-000  MV/m,  E^f  =  2.5  MV/m 


Ej~)(j  —  1.000  MV/m,  Erf  —  3.1  MV/m  Ejjq  —  1.000  MV/m,  Erf  —  2.9  MV/m 


Figure  14.  Particle  growth  for  points  near  the  lower  multipactor  boundary  curve.  In 
this  case,  /  =  2.45  GHz.  Each  of  these  points  is  shown  on  Figure  13 
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3.2.2  Space  Charge  Case. 


In  the  previous  case,  space-charge  effects  were  ignored.  While  this  is  a  good  guide 
to  what  may  happen  in  an  actual  system,  including  the  contributions  of  electrons  to 
the  restoring  held  during  the  multipactor  process  provides  for  a  more  realistic  simu¬ 
lation.  In  this  case,  Edc  is  not  constant,  but  changes  both  temporally  and  spatially, 
since  the  electron  density  also  varies  in  time  and  space.  Valfclls,  et  al.  [16]  ,  examine 
the  effects  of  the  spatial  variation  by  using  a  time-averaged  model.  Through  this 
model,  they  conclude  that  the  percentage  of  input  power  deposited  on  the  surface 
for  the  case  when  space-charge  effects  are  included  is  no  longer  a  fixed  percentage  of 
the  input  power  (as  was  determined  to  be  the  case  by  Ang,  et  al.  in  [1]),  but  the 
percentage  of  power  deposited  varies  based  on  the  magnitude  of  Erf  ■  This  result  is 
confirmed  by  Fichtl  [7]. 

In  comparing  our  simulations  with  theory  in  the  case  where  space-charge  effects 
are  taken  into  account,  it  is  more  useful  to  compare  with  the  results  presented  by 
Valfclls,  et  ah,  [16]  and  Kim  and  Verboncoeur  [9]  than  with  those  presented  by  Fichtl 
[7],  since  Fichtl  creates  susceptibility  curves  based  on  time-averaged  values  of  Edc  ■ 
While,  again,  this  is  a  useful  guide  as  to  when  multipactor  occurs  in  a  system,  it  does 
not  fully  represent  the  actual  physics  of  the  multipactor  phenomenon. 

First,  it  is  useful  to  look  at  the  time-variant  physics  of  the  problem  as  presented 
in  Verboncoeur,  et  al.  [9].  As  is  shown  in  Figure  15,  the  restoring  held,  Edc  >  varies 
periodically  in  time,  as  does  the  number  of  particles  in  the  system.  As  was  alluded 
to  in  Figure  7,  the  plot  of  EDc  versus  ERF  is  no  longer  a  straight  line.  Instead,  the 
result  is  a  parametric  plot  of  EDc  (t)  and  ERF  (t),  describing  the  system  trajectory 
in  susceptibility  space  through  the  simulation.  This  is  referred  to  as  a  ’’bowtie”  plot. 
Thus,  as  the  system  settles  into  a  periodic  steady  state,  the  values  of  ERF  and  EDc  are 
shown  to  move  along  the  periodic  curves  shown  in  Figure  16.  The  dotted  lines  are 
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created  by  the  ’’simple”  theory  (the  case  in  which  space-charge  effects  are  ignored) 
and  are  reflected  around  the  Erf  =  0  line  to  show  the  positive  and  negative  values  of 
Erf  ,  rather  than  just  the  absolute  value.  When  the  trace  is  between  the  two  upper 
and  two  lower  dotted  lines,  electron  growth  due  to  multipactor  occurs;  however,  when 
it  it  is  between  the  two  middle  dotted  lines,  the  number  of  electrons  in  the  system 
decays.  Since  the  ’’bowtie”  plot  shows  the  system  trajectory  moving  periodically 
through  regions  of  electron  density  growth  and  decay,  it  can  be  reasonably  expected 
from  these  plots,  therefore,  that  the  electron  density  changes  in  a  periodic  manner. 
The  plots  in  Figure  16  were  generated  using  the  same  parameters  given  by  Valfclls, 
which  are  listed  in  Table  1.  Consequently,  these  plots  also  confirm  Barker’s  assertion 
in  [3]  that  the  lower  boundary,  not  the  upper  boundary,  is  associated  with  saturation 
in  the  system,  since  the  traces  circulate  around  the  lower  boundaries.  This  further 
justifies  consideration  of  only  the  lower  boundary. 

Valfclls  specifies  a  seed  current  of  Jinit  =  1.3  kA/rn2,  however  he  concludes  that 
the  seed  current  is  in  fact  irrelevant  in  the  final,  steady  state  of  the  system.  This  is 
true,  for  the  most  part,  unless  the  critical  electron  density  described  in  Section  4.1.2 
is  reached,  in  which  case  total  reflection  of  the  input  power  occurs.  Additionally, 
since  the  times  for  which  dumps  are  made  are  considerably  longer  than  the  seeding 
time,  the  seed  current  does  not  greatly  affect  the  results  of  the  simulations.  In  fact, 
the  effect  of  increasing  the  seed  current  is  equivalent  to  increasing  the  macroparticle 

Table  1.  Parameters  from  Kim  and  Verboncoeur  [9]  used  to  verify  ICEPIC  simulation 
setup 


OmaxO 

2 

£max  0 

400  eV 

ks 

1 

T 

2  eV 

46 


(b)  Particles 


Figure  15.  Evolution  in  time  of  fields  and  particles  in  a  multipactoring  system,  with 
Erf  =  3  MV/m  and  /  =  1  GHz.  Note  that  the  fields  are  only  shown  for  a  small  portion 
of  the  entire  simulation,  but  the  particle  count  represents  the  simulation  in  its  entirety. 


weighting  and  decreasing  the  total  number  of  macroparticles  in  the  system.  Therefore, 
we  use  a  considerably  larger  seed  current,  in  order  to  shorten  the  time  needed  to  run 
the  simulations. 

Figure  16  confirms  the  simulations  and  theory  published  by  Kim  and  Verboncoeur 
[9].  The  shape  of  these  diagrams  correspond  very  well  with  those  given  in  those  re¬ 
sults,  however,  they  do  not  align  exactly  with  Kim  and  Verboncoeur’s  simulations. 
This  is  most  likely  due  to  three  factors:  lack  of  complete  knowledge  of  initial  param¬ 
eters  (such  as  particle  weighting  and  system  geometry);  inclusion  of  scattering  and 
reflection  effects  in  the  most  recent  simulations;  and  fundamental  differences  in  the 
ICEPIC  code  (since  it  self-consistently  calculates  all  the  Electro-Magnetic  (EM)  fields 
in  the  system).  All  of  the  simulations  with  space-charge  effects  included  demonstrated 
similar  behavior,  thereby  confirming  the  validity  of  our  input  hie. 
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(a)  Erf  =  3  MV/m,  /  =  1  GHz 


(b)  Erf  =  0.3  MV/m,  /  =  1  GHz 


(c)  =  3  MV/m,  /  =  10  GHz 

Figure  16.  Results  of  the  simulations  detailed  by  Kim  and  Verboncoeur  [9].  These 
curves  detail  the  evolution  of  Erf  and  Er>c  in  time,  representing  a  periodic  steady 
state.  The  arrows  show  the  direction  of  the  system  trajectory  and  the  areas  of  electron 
density  growth  and  decay  are  labeled. 
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3.2.3  Power  Deposition  on  the  Dielectric  Surface. 

For  the  NSC  case,  it  is  determined  by  Ang  [1],  that  the  percentage  of  power  de¬ 
posited  on  a  dielectric  surface  is  a  fixed  fraction  of  the  input  power  (~  1%).  However, 
as  is  shown  by  Valefells,  et  al.  [16],  this  is  not  the  case  when  space  charge  effects  are 
considered.  According  to  Fichtl,  in  [7],  although  ’’the  power  deposited  is  no  longer 
on  the  order  of  1  %  as  is  initially  estimated,  we  find  that  it  is  generally  less  than 
5  %  regardless  of  the  input  parameters.”  This  section  will  examine  three  different 
methods  in  which  attempts  were  made  to  verify  Fichtl’s  numerical  results  and  a  table 
of  comparative  values  is  given  at  the  end. 

3.2.3. 1  Psurf  Using  the  POWER  Dump  Files. 

It  is  not  clear  whether  Fichtl  cites  an  instantaneous  power  or  a  time-averaged 
power  when  presenting  his  results.  Since  it  is  impossible  to  directly  measure  the 
amount  of  power  deposited  in  the  dielectric  in  the  course  of  an  ICEPIC  simulation, 
Fichtl  measures  the  input  power,  the  power  transmitted  through  the  dielectric  on  the 
right  hand  side  of  the  system  and  the  power  reflected  back  through  the  left  hand  side 
of  the  system.  Then,  he  uses  the  following  formula  to  calculate  the  power  deposited 
on  the  surface  of  the  dielectric  window  under  consideration: 

Psurf  Pinput  Preflected  Ptransmitted  (30) 

Using  this  method,  Fichtl  produces  a  plot  of  the  fraction  of  power  deposited  on 
the  dielectric  window  (shown  in  Figure  17).  In  a  similar  fashion,  the  POWER  dumps 
were  used  to  obtain  the  relevant  values  for  the  simulations  performed  through  the 
course  of  this  research. 
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Fraction  of  Input  Power  Deposited  vs.  ERF 


W3 


Figure  17.  Fraction  of  input  power  deposited  on  the  dielectric  window  as  calculated  by 
Fichtl  [7] 


The  theoretical  input  power  can  be  found  by  using  the  Poynting  vector: 


6  —  2Ce°-^RFo 


Which  can  then  be  multiplied  by  the  emitting  area  to  give  the  average  input  power: 


Pin,hy  ~  Sx‘vy-v  -  2«o  ElF„x,„,yn, 


(31) 


This  theoretical  input  power  is  used  as  a  benchmark  to  ensure  the  calculations  in 
the  following  sections  are  reasonable. 


3. 2. 3. 2  Psurf  Using  PART  Dump  Files. 

While  Fichtl’s  method  is  a  perfectly  acceptable  one,  independent  validation  of 
results  is  always  desirable  for  consistency  and  additional  insight.  1CEPIC  affords  a 
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method  to  dump  the  properties  of  every  particle  that  is  collected  at  any  boundary 
in  the  system.  Due  to  the  relatively  long  length  of  the  system  in  the  ^-direction  and 
the  symmetry  in  the  x-direction,  very  few  particles  are  collected  by  the  dielectric  on 
the  left  side  of  the  system.  Therefore,  it  can  be  assumed  that  every  particle  that 
would  be  reported  by  this  dump  is  collected  at  the  boundary  of  the  dielectric  under 
study.  The  velocity  information  from  this  PART  dump  hie  can  be  used  in  conjunction 
with  the  macroparticle  weighting  information  to  calculate  the  energy  and,  in  effect, 
the  power  deposited  to  the  surface  over  a  given  dump  interval.  The  particle  dump 
hies  are  written  every  certain  number  of  time  steps,  defined  by  the  dump.int  variable 
in  the  input  hie.  In  our  case,  dump_int  was  set  to  dump  every  l/50th  of  a  period, 
where  a  period  is  17,141  time  steps  and  a  time  step  (dt)  is  ~  2.38  x  10-14  s.  The 
dump  hie  format  includes  the  position  vectors  (x,  y,  z)  as  well  as  the  velocity  vectors 
(vx,vy,vz)  for  each  of  the  collected  particles.  All  that  remains  to  be  found  at  this 
point  is  the  macroparticle  weighting,  which  is  also  included  in  the  dump  hie,  since 
the  charge  for  each  macroparticle  is  also  included  in  the  dump  hie.  Therefore,  the 


mass  of  the  macroparticle,  Mmacro,  is  found  by  using  the  macroparticle  to  electron 

ratio  ^Qmacro / Qelectroii)  • 


Qelectron 


TYl electron 


Then,  the  Energy  can  be  calculated  for  each  (ith)  particle  during  a  given  dump 
interval,  using: 


£macro,i  0  M macrovtotal,i 


where, 

2  2,2,2 

Vtotal,i  =  Vx,i  +  Vyii  +  VZii 


Since  each  dump  interval  corresponds  to  a  hxed  time  period,  the  individual  power 
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can  then  be  calculated  from: 


p ; 


macro, i 


s, 


macro,! 


dt 


\Wmev  ltali 
dt 


where  W 


Qmacro 


Qelectron 


Then  the  total  power  is  summed  for  all  the  particles  (over  i)  in  the  dump  hie,  leading 
to  a  total  power  deposited  on  the  surface  by  the  macroparticles  (Ptotai, macro)-  Since 
these  quantities  have  already  been  adjusted  to  account  for  the  macroparticle  weight¬ 
ing,  this  is  the  total  power  deposited  on  the  surface  of  the  dielectric.  Thus,  calculating 
the  power  deposited  on  the  surface  of  the  dielectric  is  a  simple  matter  of  using  the 
velocities  of  the  impinging  electrons  and  the  macroparticle  weighting  to  calculate  the 
energy  per  time  step  deposited  on  the  surface  of  the  dielectric.  Since  the  impinging 
electrons  are  the  only  source  of  power  deposited  on  the  surface,  this  appears  to  be 
the  most  accurate  and  straightforward  method  of  determining  the  power  deposited 
on  the  surface  of  the  dielectric. 


3. 2. 3. 3  Theoretical  Value  of  Psurf  Using  Valfells’  Theory. 

A  third  possible  method  of  calculating  the  power  deposited  on  the  surface  of  the 
dielectric  provides  a  theoretical  estimate.  This  method  can  then  be  compared  with 
the  previous  results  for  further  clarity  and  insight.  Upon  reexamination  of  Equations 
19  and  21,  it  is  evident  that  the  only  unknown  parameter  is  no-  However,  it  is  possible 
to  estimate  an  average  total  number  of  particles  from  the  ICEPIC  simulations  and 
solve  Equation  16  numerically,  to  obtain  a  value  of  no  (the  average  number  density 
at  x  —  0).  Care  must  be  taken  to  convert  from  macroparticles  to  real  particles  using 
the  macroparticle  to  electron  ratio  when  performing  this  calculation.  In  order  to 
calculate  no,  Equation  16  is  integrated  over  all  x  (or,  in  the  case  of  our  geometry,  z) 
and  set  equal  to  the  average  number  of  real  particles  determined  from  the  ICEPIC 
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simulation,  as  is  shown  in  Equation  32: 


W'total 

dxdy 


n  (z)  dz  =  n0 


(32) 


where  ntotai  has  been  divided  by  dxdy  in  order  to  simplify  what  is  really  a  triple 
integral  lurking  in  Equation  16  (since  the  equation  describes  a  number  density).  Then, 
Equation  32  can  be  integrated  numerically  and  no  obtained  through  a  back  solve 
process.  This  is  necessary  since  up  also  depends  on  n0.  Once  n0  is  obtained,  it  may 
be  used,  along  with  the  other  parameters  (which  are  known  and  precisely  define  the 
particular  simulation)  to  calculate  the  approximate  theoretical  power  deposition  on 
the  dielectric  surface.  As  this  is  a  mathematically  intensive  operation  which  is  to  be 
performed  over  many  iterations,  Mathematica®  software  was  used  as  an  aid.  The  full 
notebook  is  contained  in  Appendix  E.  Due  to  the  approximate  manner  in  which  the 
total  number  of  macroparticles  in  the  system  is  found,  the  result  is  relied  upon  as  a 
pure  approximation. 


3. 2. 3. 4  Comparison  of  Results  for  Power  Deposition  on  the  Dielec¬ 
tric  Surface. 

This  section  presents  a  comprehensive  comparison  of  the  three  previously  discussed 
methods  and  Fichtl’s  results.  Since  the  electron  density  varies  temporally,  it  is  clear 
that  the  power  deposited  on  the  surface  also  varies  in  time.  Therefore,  any  values 
shown  in  the  table  below  which  are  obtained  from  an  ICEPIC  simulation  are  time- 
averaged  over  a  few  periods  in  the  early  evolution  of  multipactor.  The  values  for 
Fichtl  were  obtained  by  direct  inspection  of  Figure  17.  Each  column  header  denotes 
the  specific  method  by  which  the  Psurf  (as  a  percentage  of  the  input  power)  was 
calculated.  The  final  column  is  the  no  value  in  real  particles/m3. 
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Table  2.  Comparison  of  Fichtl’s  values  for  power  deposited  on  a  dielectric  with  the 


values  found  using  the  methods  detailed  in  this  section. 


P surf 

Erf  (MV/m) 

Input  Power  (W) 

Fichtl 

POWER  dump  files 

PART  dump  files 

1.5 

9.14xl02 

0.85% 

1.78% 

2.18% 

3.0 

3.66xl03 

1.90% 

6.06% 

7.27% 

5.5 

1.23X104 

2.5% 

24.95% 

29.09% 

7.0 

1.99X104 

2.65% 

49.27% 

53.21% 

P surf 

Erf  (MV/m) 

Valfells’  theory 

no  (real  particles/m3) 

1.5 

1.50% 

4.80X1019 

3.0 

1.92% 

5.29xl02U 

5.5 

2.64% 

5.28X1022 

7.0 

2.83% 

2.12X1023 

3.3  Mitigation  Methods 

3.3.1  Thin  Coatings  of  Low-SEY  Materials. 

Adjusting  the  secondary  emission  properties  is  a  fairly  simple  task  in  an  ICEPIC 
input  file.  While  the  secondary  emission  properties  of  a  material  are  difficult  to 
measure,  Montero,  et  al.  [13]  provide  a  good  recent  reference  for  the  properties  of  thin 
TiN  films.  Some  of  the  parameters  remain  unknown  and  must  be  estimated  based 
on  an  educated  guess.  For  the  purposes  of  this  research,  the  following  parameters 
were  considered  in  simulating  a  system  in  which  the  dielectric  window  is  coated  with 
a  low-SEY  material,  such  as  TiN:0: 

The  results  of  these  simulations  will  be  discussed  in  Chapter  IV. 

3.3.2  Geometrical  Modification  of  the  Waveguide/ Window  Interface. 

Valfclls,  et  al.  [15]  investigate  the  effects  of  an  obliquely  incident  RF  field,  as 
shown  in  Figure  10.  Following  their  theoretical  development,  MATLAB®code  was 
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Table  3.  Parameters  used  to  simulate  a  thinly  coated  dielectric  surface. 


OmaxO 

1.1 

£ maxO 

300  eV 

ks 

1 

T 

2.1  eV 

developed  to  numerically  calculate  the  susceptibility  curve  for  theoretical  system  for 
a  given  angle  of  incidence,  if).  The  methodology  is  described  below,  while  the  full 
code  is  contained  in  Appendix  C.  Also,  a  method  to  simulate  oblique  incidence  using 
ICEPIC  without  modifying  the  geometry  was  developed  and  is  described  in  section 
3. 3. 2. 2. 


3.3. 2.1  Theoretical  Model  -  Matlab®for  the  NSC  Case. 

The  initialization  parameters  for  the  code  are  the  angle  of  incidence  (?/?),  vectors  of 
values  for  E^c  and  Erf  and  the  secondary  emission  parameters  (/,  £maxo,  dmaxo,  ks,  (, 
(f)  and  Vo).  The  Eoc  and  Erf  vectors  can  be  modified  to  adjust  the  ’’resolution”  of  the 
caclulations.  Note  that  E^c  and  Erf  are  scaled  to  the  dimensionless  normalization 
constants  used  by  Valfclls,  et  al.  [15].  This  is  described  in  detail  in  the  header  notes 
of  the  MATLAB®£le  contained  in  Appendix  C,  but  is  repeated  here  for  convenience: 


t 

v 

£ 

x 


treal^ 

'Ureal 

£real 


•Ureal 


yj  £max  /  nr 
Q 

TYILV  -\J  £max  /  Ufl 
U 

\/  £rna,x /uYl 


where  t  represents  time,  u  is  the  angular  frequency,  v  is  the  velocity,  £max  is  £max o  in 
units  of  Joules,  m  is  the  mass  of  an  electron  and  x  represents  the  direction  normal 
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to  the  surface.  First,  an  angle  of  incidence  is  selected.  Then,  within  an  ” inner  loop”, 
a  value  is  selected  for  Edc  ■  This  ”  inner  loop”  is  eventually  run  over  all  the  values 
contained  in  the  Edc  vector.  Then,  in  another  loop  outside  of  that  one,  the  value 
for  Erf  is  varied  in  the  same  manner.  For  each  combination  of  Edc  and  Erf  ,  the 
hoptime  is  estimated,  using  the  two  dimensional  equations  of  motion  similar  to  those 
given  in  Equation  2.2,  except  that  they  are  modified  to  account  for  the  angle  of 
incidence.  This  is  done  for  particles  emitted  at  50  different  phases  within  the  period 
of  Erf  ,  in  order  to  account  for  random  phase  emission  of  particles,  and  averaged 
over  those  phases.  In  the  following  development,  ax  is  the  ^-directed  acceleration,  ay 
is  the  ^/-directed  acceleration,  0  is  the  injection  angle  of  the  particles,  9  is  the  phase 
of  the  electric  field  at  the  time  of  emission  and  0  is  the  oblique  angle  of  incidence. 
Using  these  definitions,  the  equations  of  motion  are: 

~ax  =  Ex  =  Edc  +  Erf  sin  (t  +  9)  sin  (0) 

—ay  =  Ey  =  Erf  sin  (t  +  9)  cos  (0) 

Then,  using  the  initial  conditions,  where  vx0  and  vy0  are  the  x-  and  ^-directed  com¬ 
ponents  of  the  emission  velocity: 


Vx  (0) 

=  vx0  sin  (0) 

and 

x  (0)  =  0 

Vy  (0) 

=  vx0  cos  (0) 

and 

2/(0)  =  0 

the  complete  set  of  equations  can  be  solved,  yielding: 

vx  =  —Edc  t  +  vx0  sin  (0)  +  Erf  (-1  +  cos  (t))  cos  (9)  sin  (^) 
Erf  sin  (t)  sin  (9)  sin  (^) 

vy  =  Vxq  cos  (0)  +  Erf  (—  cos  ( 9 )  +  cos  (t  +  9))  cos  (0) 


(33) 

(34) 
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X 


(35) 


= - — - V  tvx o  sin  (0)  +  Erf  cos  (0)  (sin  (t)  -  t )  sin  (0)  - 

Erf  sin  ( 6 )  sin  (0)  +  Erf  cos  (t)  sin  ( 9 )  sin  (0) 

V  =  cos  (0)  +  -Erf  cos  (0)  (— t  cos  (0)  —  sin  ( 9 )  sin  (t  +  9))  (36) 

Since,  at  this  point  in  all  the  loops,  all  quantities  are  known,  except  for  t,  the  hoptime 
can  be  found  by  setting  Equation  35  equal  to  zero  (defined  as  the  two  points  at  which 
the  particle  is  at  the  dielectric  surface)  and  solving  for  the  corresponding  roots.  In  this 
code,  the  hoptime  is  obtained  numerically,  using  the  fzeroO  function  in  Matlab®. 
Once  it  is  obtained,  the  hoptime  is  used  to  calculate  (,  5,  Si,  S-2  and  £max  from  the 
equations  of  motion.  The  secondary  emission  coefficient  (h)  is  averaged  over  the  50 
emission  phases  for  each  pair  of  Erf  and  E^c ,  yielding  5aVg •  If  the  value  for  5avg  is 
less  than  one,  it  is  plotted  as  a  blue  circle  and  if  it  is  greater  than  one,  it  is  plotted  as 
a  red  plus.  The  rest  of  code  is  then  looped  for  many  different  angles  of  obliqueness. 
The  results  of  these  calculations  are  given  in  Chapter  IV. 

3.3.2. 2  ICEPIC  Model. 

To  date,  no  ICEPIC  simulations  of  obliquely  incident  fields  have  been  performed 
to  verify  the  narrowing  of  the  multipactor  boundary  curves.  Because  of  the  limitations 
of  the  commands  that  must  be  used  to  create  shapes  in  ICEPIC  input  hies,  creating 
a  new  geometry  for  these  ICEPIC  simulations  would  be  incredibly  complex  and  time- 
consuming.  Therefore,  the  time  variable  capabilities  contained  within  ICEPIC  were 
applied  to  the  Er>c  held  to  simulate  the  interaction  of  the  Erf  held  with  the  E^c  held, 
in  the  NSC  case.  The  theory  behind  this  simplification  is  relatively  straightforward. 
According  to  Figure  10  the  Erf  held  is  described  by: 

Erf  =  Erfo  sin  (ut  +  9) 
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However,  in  the  case  of  an  obliquely  incident  RF  field,  Erf  takes  the  form  of: 

ErF  =  ERF0  sin  (ut  +  9)  sin  (-0) 

As  electrons  leave  the  surface,  Edc  evolves  as  a  ’’constant”  held.  Therefore,  the  total 
held  in  the  ^-direction  is  given  by  (also  assuming  the  simulation  starts  when  9  —  0): 

EZytot  =  Edc  +  Erfo  sin  (ut)  sin  ($) 

This  leads  to  a  time- varying  EDc  held,  which  is  predicted  by  Valfells,  et  al.  [15]. 
This  time  varying  Edc  will  naturally  lead  to  a  change  in  the  susceptibility  of  the 
system  to  multipactor.  It  is  worth  noting,  though,  that  previous  research  did  not 
consider  the  time-varying  aspect  of  multipactor  evolution,  but  dealt  purely  in  aver¬ 
ages.  In  our  simulations,  the  time-varying  aspect  is  unavoidable  and  so  the  results 
should  be  compared  with  caution.  A  sample  input  hie  is  given  in  Appendix  D,  where 
the  time  variable  parameters  Co  and  C\  are  used  to  create  the  time  varying  Edc  held. 
According  to  the  ICEPIC  manual,  the  time  variables  can  be  used  to  create  a  held 
described  by: 

Edc, =  Co  +  Ci  sin  (ut  +  9)  +  ... 

where: 

C0  =  1  and  Ci  =  (ERF o/  -  EDC  )  sin  (^) 

Figure  18  compares  the  theoretical  output  with  actual  output  from  a  corresponding 
ICEPIC  simulation.  The  plot  in  the  theoretical  case  has  been  shifted  to  the  right  to 
account  for  the  propagation  delay  in  the  simulation.  This  is  due  to  the  fact  that  the 
wave  is  launched  from  the  left  hand  dielectric,  propagates  for  a  finite  amount  of  time 
and  is  measured  at  the  right  hand  dielectric  surface.  Also,  in  the  plots  of  the  helds 
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from  the  simulation,  a  ’’ramp  up”  time  is  noticeable  for  both  the  Erf  and  E/x?  fields 
to  move  to  the  expected  values,  which  is  a  function  of  the  time  variable  configuration 
in  ICEPIC  and  unavoidable. 


Erf  &  Edc  ( Erf  =  1.5  [MV/m],  EDC  =  0.125  [MV/m], 


=  5°) 

Erf.  ICEPIC 

- Edc ,  ICEPIC 

- Erf,  theory 

Edc ,  theory 


Figure  18.  Comparison  of  the  theoretical  and  ICEPIC  Edc  and  Erf  fields  in  the  oblique 
case.  The  theoretical  case  has  been  run  out  to  longer  times  than 
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IV.  Analysis  of  Results 


This  chapter  looks  at  the  results  of  the  results  from  the  calculations  and  simula¬ 
tions  described  in  Chapters  If  -  111  and  presents  some  analysis  regarding  these  results, 
relating  them  to  the  comprehensive  purpose  of  multipactor  mitigation  in  HPM  sys¬ 
tems.  The  destructive  mechanism  of  multipactor  is  twofold.  First  and  foremost,  the 
continual  bombardment  of  the  dielectric  surface  by  high  energy  electrons  results  in 
the  deposition  of  a  percentage  of  the  input  power  on  the  surface  of  the  dielectric. 
This  power  deposition  eventually  results  in  damage  to  the  dielectric  surface,  even¬ 
tually  causing  system  failure.  This  limits  the  input  power  of  an  HPM  system.  Also 
of  interest  is  how  the  field  strength  of  this  limit  compares  with  the  breakdown  field 
strength  of  air,  which  is  a  fixed  threshold.  Secondly,  the  high  density  of  electrons 
reflects  a  portion  of  the  Erf  field  being  reflected  back  towards  the  source  in  a  real 
HPM  system.  Eventually,  this  can  results  in  damage  to  or  failure  of  the  source. 
Both  of  these  mechanisms  are  discussed  below,  followed  by  the  effect  of  the  proposed 
mitigation  strategies. 

4.1  Multipactor  Power 

The  power  absorbed  by  the  electrons  of  a  multipactoring  system  can  be  viewed  in 
many  different  ways.  As  was  previously  stated,  the  power  deposition  on  the  surface 
of  the  dielectric  is  of  primary  interest.  Also  of  concern  is  the  power  reflected  back 
towards  the  source  by  the  evolving  electron  density,  which  can  result  in  damage  to 
the  HPM  source.  Also,  whether  or  not  multipactor  supercedes  the  limit  imposed  by 
air  breakdown  is  considered  in  this  section. 
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4.1.1  Power  Deposition. 


As  is  evident  from  Table  2,  the  PART  dump  and  POWER  dump  methods  of  deter¬ 
mining  power  deposited  on  the  surface  of  the  dielectric  correlated  well,  but  produced 
drastically  different  results  than  the  method  using  Valfells’  theoretical  method  and 
Fichtl’s  results.  This  remains  an  unreconciled  discrepancy.  The  two  different  dump 
files  used  from  ICEPIC  simulations  are  independent  of  one  another  and  are  well 
tested.  Also,  the  input  power  from  the  dump  hies  match  up  almost  exactly  with  the 
theoretical  input  power,  as  calculated  in  Equation  31.  Since  the  PART  dump  method 
uses  basic  physics  to  calculate  the  kinetic  energy  and  power  associated  with  a  particle 
impacting  the  dielectric  surface  over  a  short  time,  this  seems  to  be  the  simplest  and 
most  trustworthy  method.  However,  it  is  possible  that  ICEPIC  contains  a  code  bug 
that  incorrectly  reports  values  to  either  both  of  these  dump  hies,  or  a  power  sink 
exists  that  has  been  overlooked.  Fichtl  [7],  considers  later  times,  but  this  would  rep¬ 
resent  an  increase  in  the  amount  of  power  deposited  on  the  dielectric  surface,  since 
the  electron  density  would  also  be  higher.  It  is  unclear  what  these  results  signify  and 
how  they  should  be  interpreted.  Further  study  is  recommended  to  correlate  these 
results  to  the  literature. 

4.1.2  Reflected  Power. 

In  the  course  of  performing  ICEPIC  simulations,  an  interesting  phenomenon  came 
to  light.  For  many  simulations,  at  lower  values  of  ErtF  ,  the  amount  of  reflected  power 
was  on  the  order  of  tenths  of  a  percent  of  the  input  power.  However,  as  both  the 
particle  weighting  and  value  of  ERF  were  increased,  and  the  simulations  were  allowed 
to  run  out  to  longer  times,  the  reflected  power  was  seen  to  greatly  increase.  While 
increasing  the  particle  weight  was  geared  to  decrease  the  simulation  time,  it  had  the 
unintended  result  of  producing  very  high  real  electron  densities  at  these  later  times. 
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In  fact,  at  electron  densities  on  the  order  of  5  x  1019  electrons/m3  ,  a  near  complete 
reflection  of  input  power  occurred.  This  is  not  accounted  for  in  any  of  the  current 
theories.  It  is  possible  that  this  phenomenon  could  account  for  the  ’’knee”  Fichtl 
notes  in  his  power  deposition  analysis.  The  power  deposited  on  the  surface  should 
rise  proportional  to  the  Erf  ,  however,  when  an  increasingly  larger  amount  of  power  is 
reflected,  the  power  deposited  on  the  surface  would  decrease  accordingly.  Obviously, 
reflected  power  is  of  great  concern  to  preserving  expensive  and  sensitive  HPM  sources. 

4.1.3  Multipactor  Power  Limit  versus  Air  Breakdown  Limit. 

Ambient  air  breakdown  is  a  fixed  limit,  one  that  truly  restricts  every  HPM  system. 
In  the  case  where  SC  is  considered,  the  susceptibility  curve  is  no  longer  interpreted  as 
a  metric  of  whether  or  not  multipactor  occurs,  but  it  describes  the  system  trajectory 
in  susceptibility  space  as  it  moves  through  regions  of  electron  growth  and  decay. 
However,  there  are  still  some  values  of  Erf  for  which  multipactor  does  not  occur, 
because  the  RF  held  does  not  impart  enough  energy  to  the  electrons  to  facilitate 
multipactor  growth.  Where  this  boundary  falls  in  relation  to  the  fixed  limit  of  ambient 
air  breakdown  holds  drastic  implications  for  the  future  of  higher  powered  systems.  A 
series  of  ICEPIC  simulations  was  performed  to  find  this  lower  boundary  and  the  value 
of  Erf  at  which  a  system  no  longer  sustains  multipactor  growth  is  approximately  0.25 
MV/m.  The  held  breakdown  strength  of  ambient  air  at  sea  level,  for  a  reasonably 
assumed  emitting  area  of  1  m2  is  given  to  be  2.745  MV/m  [2],  In  the  context  of 
the  ’’bowtie”  plots  (see  Figure  16),  the  behavior  correlates  well  with  our  current 
understanding  of  multipactor.  As  the  overall  amplitude  of  ERF  decreases,  the  total 
number  of  particles  in  the  system  decreases  as  well.  Therefore,  both  the  height  and 
width  of  the  ’’bowtie”  decrease,  to  the  point  where  the  ’’bowtie”  no  longer  makes 
excursions  into  the  region  of  multipactor  growth  and  all  the  particles  in  the  system 
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are  lost  to  the  dielectric  surface.  Therefore,  since  the  multipactor  limit  is  lower  than 
the  ambient  breakdown  limit,  it  is  evident  that  this  is  the  reason  why  multipactor 
remains  a  serious  concern  in  the  design  and  implementation  of  HPM  systems.  The 
upper  limit  was  found  to  be  greater  than  1  x  109  V/m,  which  is  well  above  the 
capabilities  of  modern  systems,  so  the  upper  limit  is  not  considered. 

4.2  Mitigation  Techniques 

4.2.1  Geometric  Modification  of  the  Waveguide/ Window  Interface. 

Using  the  code  described  in  Section  3. 3. 2.1,  susceptibility  curves  for  a  certain 
system  were  generated  for  the  NSC  case.  These  curves  are  shown  in  Figure  19.  From 
these  curves,  it  is  found  for  smaller  angles  of  incidence,  0  ~  5  — 10°,  the  susceptibility 
region  actually  increases.  However,  for  angles  of  incidence  greater  than  or  equal  to 
10°,  the  susceptibility  region  begins  to  dissipate  and  then  contracts  rapidly.  At  angles 
greater  than  20°,  the  system  appears  to  be  no  longer  susceptible  to  multipactor  at  all. 
It  is  not  clear  from  this  research  whether  or  not  the  susceptibility  actually  disappears 
or  if  this  is  numerical  artifact.  Regardless,  the  diminished  region  of  susceptibility  is 
promising  and  confirms  the  theory. 

With  regards  to  the  ICEPIC  simulations,  the  results  are  not  as  conclusive  as  with 
the  MATLAB®calculations.  The  results  of  the  simulations  are  compared  with  the 
MATLAB® results  in  Figure  20.  From  these  simulations,  it  is  clear  that  the  results 
correlate  well  with  the  theory  and,  in  fact,  show  the  lower  boundary  moving  up.  It  is 
unclear  what  the  SC  case  might  produce  in  ICEPIC  simulations.  One  can  imagine, 
though,  from  looking  at  Figure  16,  that  since  the  susceptibility  range  is  significantly 
reduced,  the  system  trajectory  may  spend  less  time  in  the  range  which  represents 
multipactor  growth  and,  therefore,  reach  a  periodic  steady  state  at  a  lower  average 
electron  density.  Therefore,  one  might  anticipate  that  less  power  would  be  deposited 
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on  the  dielectric  surface  and  less  power  would  be  reflected  towards  the  source.  This, 
in  turn,  would  allow  systems  to  operate  at  higher  powers  before  reaching  the  point 
where  enough  power  is  deposited  on  the  surface  to  cause  failure  of  the  dielectric. 

4.2.2  Thin  Coatings  of  Low-SEY  Materials. 

Reducing  the  secondary  emission  coefficient  of  the  dielectric  surfaces  has  the  pri¬ 
mary  effect  of  reducing  the  electron  density  in  the  system,  which  gives  a  twofold 
advantage:  a  reduction  in  the  power  deposited  on  the  surface  of  the  dielectric  and 
a  reduction  in  the  electron  density  of  the  system  which  equals  less  power  reflected 
towards  the  source.  Figure  21  compares  the  power  deposited  on  the  surface  of  the 
dielectric  for  the  same  four  cases  run  in  Table  2.  It  is  evident  from  Figure  21  that 
the  power  deposited  on  the  surface  of  a  dielectric  is  drastically  reduced  for  a  system 
coated  in  a  him  like  Ti:N. 
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Figure  19.  Susceptibility  curves  generated  from  Mat  lab®  simulations,  where  Edc  and 
Erf  are  scaled  as  [MV/m]  x  (//1GHz)-1  (£maa;o/400eV)  1^2.  In  this  system,  Smax o  =  3. 
Plus  signs  (+)  represent  points  at  which  multipactor  occurs  and  circles  (o)  represent 
areas  where  multipactor  does  not  occur. 
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(a)  ip  =  0° 


(b)  ip  =  10° 


Edc  ,  scaled 

(c)  ip  =  15° 


Figure  20.  Susceptibility  curves  generated  from  ICEPIC  simulations,  where  Eqc  and 
Erf  are  scaled  as  [MV/m]  x  (//  1GHz)  1  (fmaa:o/400eV)_1^2.  In  this  case,  only  the  lower 
boundary  is  shown. 
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Fraction  of  Input  Power  Deposited  on  Surface 


Figure  21.  Comparison  of  the  power  deposition  metrics  for  the  bare  dielectric  and  the 
dielectric  covered  by  a  thin  coating  of  a  low-SEY  material. 
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V.  Conclusions  and  Future  Work 


5.1  Conclusions 

Chapter  II  presented  the  basic  theory  of  multipactor  and  offered  a  short  summary 
of  the  known  theory  relating  to  this  phenomenon.  Secondary  Electron  Emission  (SEE) 
in  materials  is  the  underlying  phenomenon  leading  to  a  net  gain  in  electrons,  or  a 
SEY  greater  than  one  and  an  eventual  electron  avalanche  in  HPM  systems.  Secondary 
emission  is  well  described  by  Vaughan  in  [17]  and  is  dependent  on  a  few  parameters 
of  the  impacted  material,  as  well  as  the  energy  of  the  impacting  electrons.  Both 
the  NSC  and  SC  cases  were  disucussed,  along  with  the  different  interpretations  of 
each  case.  In  the  NSC  case,  ’’universal”  multipactor  boundaries  may  be  constructed, 
normalized  to  Emax0  and  /,  so  that  multipactoring  regimes  may  be  predicted  for 
many  systems.  For  the  SC  case,  the  electron  density  and,  accordingly,  E Dq  varies 
both  spatially  and  temporally,  so  the  theory  differs  and  the  susceptibility  curves 
must  be  also  be  interpreted  according  to  the  temporal  variation.  In  light  of  this 
discovery,  the  term  ”Ejjc”  is  misleading,  since,  in  reality,  the  restoring  held  is  not 
constant.  However,  in  order  to  be  consistent  with  the  literature  and  current  jargon, 
the  term  E^c  is  still  used.  So-called  ’’bowtie”  diagrams  can  be  created  from  the 
relationships  between  Erf  and  Eoc  and  the  guidlines  from  the  NSC  case  estimate 
when  the  secondary  emission  coefficient  is  greater  than  one  and  when  it  is  less  than 
one.  An  overview  of  general  PIC  code  theory  and  input  parameters  served  as  the 
foundation  for  understanding  ICEPIC  results  and  how  secondary  emission  effects  are 
included  in  ICEPIC  simulations.  With  regard  to  mitigation  methods,  two  methods 
were  discussed  in  this  paper  -  thin  coatings  of  low-SEY  materials  and  RF  held  oblique 
incidence.  The  theory  behind  each  was  discussed. 

In  Chapter  III,  the  rationale  behind  the  selected  input  parameters  was  discussed, 


along  with  a  detailed  description  of  important  variables  in  ICEPIC.  Then,  current 
simulations  were  compared  with  those  presented  in  the  literature.  ICEPIC  simula¬ 
tions  confirmed  the  results  of  both  Verboncoeur  [9]  and  Fichtl  [7],  while  updates  to  the 
code  resulted  in  small  changes  in  the  universal  susceptibility  curves.  Current  ICEPIC 
simulations  show  that,  in  some  cases  (values  of  Erf  less  than  3.0  MV/m),  at  rela¬ 
tively  short  times  (~  5  periods),  the  particle  density  reaches  a  periodic  steady  state. 
In  some  cases,  though,  the  particle  density  continues  to  rise  slowly  throughout  the 
simulation  and  may,  at  longer  times,  reach  a  similiar  steady  state.  For  the  purposes 
of  this  paper,  the  early  evolution  of  multipactor  is  studied,  where  values  for  power 
and  holds  are  reported  between  five  and  eight  periods.  While  this  may  not,  in  every 
case,  represent  a  true  ” steady  state,”  it  is  taken  to  accurately  represent  a  snapshot  in 
time  of  the  order  of  magnitude  of  multipactor  and  can  be  extended  to  longer  times. 
Power  deposition  on  the  surface  is  perhaps  the  most  destructive  effect  of  multipactor, 
resulting  in  failure  of  the  dielectric  window.  Therefore,  a  large  amount  of  time  was 
spent  understanding  and  analyzing  power  deposition.  Three  different  methods  were 
used  to  estimate  the  power  deposited  on  the  surface.  Two  of  the  methods,  inde¬ 
pendent  dumps  from  the  ICEPIC  simulations,  confirmed  each  other,  while  the  third, 
implementing  the  space  charge  theory  from  [16]  correlated  well  with  Fichtl’s  results 
in  [7].  The  reason  for  the  large  disparity  between  the  two  different  results  remains 
unreconciled.  Since  Fichtl’s  input  hies  were  the  basis  for  the  simulations  presented  in 
this  paper  and  all  known,  realistic  power  sinks  were  accounted  for,  no  determination 
can  be  made  yet  as  to  the  source  of  this  disparity.  Additionally,  the  ICEPIC  models 
behind  the  two  mitigation  methods  were  presented,  and  the  computational  theory 
behind  development  of  susceptibility  curves  in  MATLAB  explained. 

It  has  been  demonstrated  that  multipactor  growth  occurs  very  quickly  in  a  system, 
usually  within  a  few  nanoseconds.  Considering  that  a  pulse  in  an  HPM  system 
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may  be  on  the  order  of  100  ns,  it  is  clear  that  a  large  amount  of  power  may  be 
deposited  on  the  surface  of  the  dielectric  in  a  very  few  pulses.  Analysis  of  the  results 
of  power  deposition  metrics  and  simulations  for  the  mitigation  methods  was  presented 
in  Chapter  IV.  It  was  found  that  the  electron  density  is  important  not  only  because 
of  the  power  deposited  on  the  surface  of  the  dielectric,  but  because  higher  densities 
reflect  large  amounts  of  power  back  towards  the  source  and  may  results  in  damage  to 
the  HPM  source. 

For  the  SC  case,  lower  and  upper  boundaries  were  founnd  -  outside  of  which  multi- 
pactor  no  longer  occurs  within  a  system.  The  breakdown  field  strength  of  ambient  air 
is  between  these  two  values,  emphasizing  the  need  to  consider  multipactor  evolution 
and,  most  certainly,  mitigation  methods,  when  designing  an  HPM  system. 

Also  presented  in  Chapter  IV  is  an  analysis  of  the  the  effects  RF  obliqueness 
and  on  the  susceptibility  curves.  The  curves  are  significantly  narrowed  for  angles  of 
obliqueness  greater  than  10°.  Since  it  is  primarily  the  upper  boundary  that  shifts 
when  considering  an  obliquely  incident  RF  field  and  it  is  the  lower  boundary  that  is 
primarily  associated  with  electron  saturation,  the  goal  would  be  to  create  a  system 
in  which  the  angle  of  obliqueness  is  greater  than  20°,  so  that  the  susceptibility  to 
multipactor  would  be  very  small.  It  was  also  found  that,  clue  to  a  narrowing  of  the 
susceptibility  region,  in  a  carefully  designed  system,  less  power  may  be  deposited  on 
the  dielectric.  Therefore,  systems  with  higher  input  powers  may  be  developed. 

In  general,  it  is  found  that  both  mitigation  methods  that  were  explored  in  this 
paper  result  in  a  dampening,  and,  for  some  systems,  total  elimination,  of  electron 
multipactor  growth  and  a  decrease  in  the  amount  of  power  deposited  on  the  surface 
of  the  dielectric.  Therefore,  a  combination  of  these  two  methods  deserves  further 
study  through  simulation  and  consideration  when  designing  current  and  future  HPM 
systems. 
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5.2  Future  Work 


Many  aspects  of  the  work  presented  in  this  paper  require  further  work  in  order 
to  justify  their  inclusion  in  modern  HPM  systems.  Some  of  these  aspects  are  natural 
extensions  of  the  research  described  in  this  paper,  while  others  are  required  to  further 
solidify  the  correlation  between  simulation  and  theory.  Suggestions  for  future  work 
are  made  with  the  goal  of  obtaining  a  more  complete  and  physically  accurate  model 
for  simulating  the  multipactor  phenomenon  and  the  possible  mitigation  techniques 
discussed  in  this  paper  (as  well  as  future  techniques). 

The  simulations  that  were  run  during  the  course  of  this  research  were  very  spe¬ 
cialized;  they  considered  only  the  multipactor  effect  inside  the  system  for  generic, 
widely-accepted  and  utilized  parameters.  However,  it  is  not  known  how  the  pa¬ 
rameters  of  current  simulations  compare  with  more  realistic  parameters  (i.e.  size, 
emissivity)  of  real  dielectric  windows.  One  possible  area  for  future  work  is  to  update 
these  parameters  and  analyze  the  effects  of  multipactor  in  a  more  realistic  sense. 

While  some  simulations  of  a  system  in  which  the  dielectric  was  coated  with  a 
thin,  low-SEY  material  were  conducted,  a  few  assumptions  were  made  in  running 
these  simulations.  Each  of  these  assumptions  should  be  studied  and  tested,  so  as 
to  be  proven  or  discarded  in  future  simulations.  Most  importantly,  it  was  assumed 
that  the  SEY  coefficient  of  the  dielectric  surface  is  reduced  to  that  of  the  thin  film 
covering.  However,  it  is  unclear  if  this  is  actually  the  case.  It  is  possible  that  some 
electrons  strike  the  surface  and  penetrate  to  the  original  dielectric.  Obviously,  this  is 
highly  dependent  on  the  thickness  of  the  thin  film.  Therefore,  the  average  secondary 
coefficient  may  be  reduced  by  some  amount  less  than  was  assumed  in  this  paper.  A 
model  could  be  developed  which  includes  the  physics  of  multiple  layers  of  materials 
with  different  SEY  coefficients.  This  model  could  then  be  used  to  study  the  necessary 
thickness  of  the  coating  of  a  low-SEY  material. 
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This  paper  describes  a  few  simulations  that  were  run  to  emulate  an  obliquely 
incident  RF  held  for  the  NSC  case  by  using  the  time  variables  inherent  in  ICEPIC. 
This  is  useful  in  determining  the  multipactor  boundaries  as  presented  by  Kishek,  et 
al.  [11],  However,  the  power  deposition  on  the  surface  cannot  be  studied  accurately 
in  this  case,  since  the  space  charge  effects  are  ignored.  Therefore,  it  would  be  useful  to 
run  these  simulations  in  the  SC  case.  This  would  require  altering  the  geometry  of  the 
system  to  include  a  window  that  includes  an  RF  window  that  is  physically  oblique  to 
the  incoming  RF  held.  This  would  be  a  complex  undertaking  in  ICEPIC,  due  to  many 
of  the  unknown  factors,  two  of  which  are:  actually  creating  the  geometry  in  ICEPIC 
and  understanding  how  the  PML  boundaries  handle  this  type  of  geometry,  since  it 
now  appears  asymmetric  to  the  PML.  Implementing  a  truly  oblique  geometry  in 
ICEPIC  would  allow  simulation  and  analysis  of  many  types  of  multipactor-inhibiting 
geometries,  to  more  accurately  determine  the  effectiveness  of  this  mitigation  strategy. 

According  to  [3],  outgassing  effects  at  the  surface  of  the  dielectric  produce  an 
as  yet  unknown  quantitative  variation  in  the  SEY  of  the  system.  Since  impacting 
electrons  will  cause  desorption  of  gas  at  the  dielectric  surface,  at  certain  levels  of 
Erf  ,  a  layer  of  gas  may  form  above  the  surface  of  the  dielectric  (a  process  in  itself 
not  well  understood)  and  potentially  form  a  plasma,  resulting  in  near  total  reflection 
of  the  incident  power.  The  effect  of  this  gas  layer  on  electron  multipactor  in  the 
system  is  not  well  studied  and  deserves  further  attention. 
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Appendix  A.  Sample  ICEPIC  Input  File  for  Non-Space 

Charge  Case 


This  input  file  is  for  the  non-space  charge  case. 

1  [Defaults] 

ENERGY=1.5e6  ;  these  parameters  can  be  varied  for  surveys 
EDC  =0  .  1 25  e6 

;  this  is  fichtl’s  original  file  case  and  is  meant  to  recreate  his... 
thesis  work  . 

6  ;  it  has  been  updated  to  work  with  current  versions  of  ICEPIC. 

[Variables] 
c0=299792458 . 0 
pi =3 . 141592653589793 
11  eps0=8 . 85418781762E-12 
me  =9 . 1093897E-31 
qe  =  1 . 60217738E-19 
test  =  1 
f c=2 . 45E9 
16  Iamda0=c0/fc 
ez=EDC 

energy = ENERGY 

dz=lamda0 /3000  ;  determines  the  z-resolution  of  the  system 

PMLdepth=int ( ( . 025* lamdaO ) / dz+ . 5) 

21  Nx=50 
Nz  =300 
xmin=0 
xmax=Nx*dz 
zmin=0 

26  zmax=Nz*dz 

dt=0.99*dz/ (c0*sqrt (2) ) 
period=int (1/fc/dt) 
max step  =  15*  period 

hopt ime  =  (2*me / qe / ez )* sqrt (2* qe *420/me ) *0 . 8509035245  ;  this  is  ... 

worst  -  case 

31  hopdist=sqrt (2*qe*420/me)*hoptime 
Ndielect=10*dz 

print  "Hop  time=",  hoptime  ; comments  in  the  * . dat  file 

print  "Hop  distance=",  hopdist 
print  " TimeSt eps /Hopt ime = " ,hoptime/dt 
36  print  "Hop  distance /dz=" ,  hopdist/dz 

[Cartesian] 

ig_dir  =  (0 ,1,0) 

def ault_e = (0 . ,0. , -ez) 

41 

[XGridl] 

range=Unif orm (xmin -3*dz , xmax+3*dz , Nx+6)  ;  accomodates  symmetry  ... 

shortcomings 
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[YGridl] 

46  range=Unif orm (0 , 1 , 1 ) 

[ZGridl] 

range =Uni form (zmin , zmax , Nz) 

51  [Symmetry] 

; symmetry  =  ( xmin , - 1 .  ,zmin-10*dz, xmax ,  2  .  , zmax+10*dz) 
periodicl=(X , xmin , xmax ) 

;  sym  mode  =  (  1 . 0  ,  0 . 0  ,  0 . 0  ) 


56 

[Time] 
dt  =  dt 

step_max=maxstep 
courant_value= . 95 
61  kill_if _below_part s =5 
kill_if_above_parts=2e6 

[ShapeN ] 

;  computational  domain 
66  shape=Box (xmin , xmin , zmin , xmax , xmax , zmax) 

[ShapeN ] 

; Dielectric  on  RHS  of  system 

shape=Box(xmin ,0.0 , zmax -4*dz-PMLdepth*dz-Ndielect ,xmax ,1.0, zmax -4* 
dz-PMLdepth*dz) 

71  material=dielectric 

[ShapeN ] 

; Dielectric  on  LHS  of  system 

shape=Box(xmin ,0.0 , zmin+4*dz+PMLdepth*dz ,xmax , 1 .0 , zmin+PMLdepth*dz 
+Ndielect+4*dz) 

76  mater ial=dielectr ic 

[MaterialN] 
name=dielectric 
epsilon=l . 001 

81 

[ParticlesN] 

; Secondary  emitter  on  dielectric  on  RHS  of  system 

shape  =  Box (xmin-dz ,0.0,zmax-4.25*dz-PMLdepth*dz-Ndielect  ,xmax  ,1.0,  . 

zmax -3 . 75*dz-PMLdepth*dz-Ndielect ) 
method = SECONDARY 
86  se c_coef f i c i ent =3 
sec_energy  =2 . 1 
sec_threshold=l 
sec.refl =0 . 03 
sec.scat =0 . 07 

91  sec_energy_max_yield=420 
sec_ks=l 
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96 


q=-qe 
mass  =me 
dir  =2 

[ParticlesN] 

;Particle  emitter  on  dielectric  on  RHS  of  system 
shape=Box(xmin-dz,0.0,zmax-4.25*dz-PMLdepth*dz-Ndielect 
zmax -3 . 75*dz-PMLdepth*dz-Ndielect ) 
method=BEAM 
101  temp =2 . 1/ test 

energy =0 . 0/test 
random=l 
q=-qe/test 
mass=me/ test 
106  current =1000/ test 

inject_interval=100 
flattop=.5/fc 
amp2  =0 
tstop= . 5/f c 
111  smooth=l 
dir  =2 

[PlanewaveN ] 

; Generates  Planewave  from  xmin  to  xmax  across  system 
116  shape=Box(xmin-2*dz,0.0,zmin+PMLdepth*dz+2*dz,xmax+2*dz 
PMLdepth*dz-2*dz) 
theta=180 . 0 
phi =0 . 0 
psi =180 . 0 
f requency=f c 

121  origin=(xmin ,0.0 , zmin+2*dz+dz*PMLdepth) 

E0=energy 

[Expert ] 

; cell_mult =12  ; not  sure  what  this  parameter  is,  not  in 

ICEPIC  version 

126  r eque st ed_max_num_part icles=3e6 

self _f i elds =0  ;this  is  what  turns  off  the  space  charge 

[BoundN ] 

;  PML  on  LHS  of  system 
131  method=PML 

depth=PMLdepth 
R=1 . 0E-4 
t  aper  =0 
order  =2 
136  dir=2 

shape=Box (xmin ,0 , zmin-dz/2 , xmax , 1 . 0 , zmin+dz/2) 

[BoundN ] 

;  PML  on  RHS  of  system 
141  method=PML 


xmax  ,1.0, 


1.0, zmax - 


current 
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depth=PMLdepth 
R=1 . OE-4 
t  aper  =0 
order  =2 
146  dir =2 

shape  =  Box ( xmin ,0 , zmax - dz /2 , xmax , 1 . 0 , ziax  +  dz/2) 

;  no  dumps  for  this  case,  since  all  we’re  looking  for  is  mp  growth 
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Appendix  B.  Sample  ICEPIC  Input  File  for  Space  Charge 

Case 

This  input  file  is  for  the  space  charge  case. 

[Def  ault  s ] 

ENERGY=3e6  ;  these 
EDC  =0 

;  this  is  fichtl’s  original  file  and  is  meant  to  recreate  his  ... 
thesis  work  in  the 

;  space  charge  case .  it  has  been  updated  to  work  with  current  . .  . 
version  of 

;  ICEPIC.  put  in  new  field  dumps  to  correspond  with  the  power  ... 
dumps  and 

;  shortened  the  max  time  of  the  simulation. 

) 

[Variables] 

c0=299792458 . 0 

pi =3 . 141592653589793 

eps0=8 . 85418781762E-12 

me  =9 . 1093897E-31 

qe  =  1 . 60217738E-19 

test  =  1 

f c=2 . 45E9 

Iamda0=c0/ f c 

ez=EDC*le6 

energy = ENERGY 

dz= lamdaO / 12000  ;  finer  resolution  than  the  NSC  case 

PMLdepth=int ( ( . 0125* lamdaO ) / dz+  .  5) 

Nx  =  3 

Nz =200+2* PMLdepth 

xmin=0 

xmax=Nx*dz 

zmin=0 

zmax  =  Nz  *  dz 

dt=0 . 99*  dz / (c0*sqrt (2) ) 
zbuf  f  =6*  dz 
period=int (1/fc/dt) 
max step  =  8*  period 
Ndielect=10*dz 
numstart  =  5* period 
numstop  =  8* period 
Dumplnt =period/50 

print  "PMLdepth=",  PMLdepth  ;  comments  in  the  * . dat  file 

print  " dz=" ,  dz 

print  " dt  =  "  ,  dt 

print  "Period=",  period 

print  "last  step=",  maxstep 
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[Cartesian] 
ig_dir  =  (0 ,1,0) 

46  def ault_e  =  (0 .  ,0.  ,-ez) 

[XGridl] 

range =Uni form (xmin -1 . 0*dz , xmax+2 ,0*dz,Nx+3) 

51  [YGridl] 

range =Uni form (0,1,1) 

[ZGridl] 

range  =  Unif  on  (  zmin  ,  zmax  ,  Nz  ) 

56 

[Symmetry] 

;  symmetry  =  (  xmin  -  .  25*  dz  ,  -1  ,  zmin  -  1 0*  dz  ,  xmax  +  .  25*  dz  ,  2  ,  zmax  + 10*  dz  ) 

;  sym  mode  =  (  1 . 0  ,  0 . 0  ,  0 . 0  ) 

periodicl  =  (X,xmin-.25*dz, xmax  + . 25*  dz ) 

61 

[Time] 
dt  =  dt 

step_max=maxstep 
courant_value= . 95 

66 

[ ShapeN ] 

shape=Box (xmin , xmin , zmin , xmax , 1 . 0 , zmax)  ;  the  computational  ... 
domain 

[ShapeN ] 

71  ; Dielectric  on  RHS  of  system 

shape=Box(xmin ,0.0 , zmax -2*zbuff -PMLdepth*dz-Ndielect , xmax , 1 . 0 , zmax 
-2*zbuff -PMLdepth*dz) 
material=dielectric 

[ShapeN ] 

76  ; Dielectric  on  LHS  of  system 

shape=Box(xmin ,0.0 , zmin+2*zbuff+PMLdepth*dz , xmax , 1 . 0 , zmin+PMLdepth 
*dz+Ndielect+2*zbuff) 
material=dielectric 

[MaterialN] 

81  name=dielectric 
epsilon=l . 001 

[ParticlesN] 

; Secondary  emitter  on  dielectric  on  RHS  of  system 
86  shape=Box (xmin -dz , 0 . 0 , zmax -2*zbuff -0 . 25*dz-PMLdepth*dz-Ndielect , . . 
xmax , 1 . 0 , zmax -2*  zbuf  f  +0 . 25*dz-PMLdepth*dz-Ndielect ) 
method = SECONDARY 
sec_coefficient=3 
sec_energy=2 . 1/ test 
sec_threshold=l/ test 
91  sec.ref 1=0 . 03 
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sec_scat  =0 . 07 

sec_energy_max_yield=420/ test 
sec_ks=l 
q=-qe/test 
96  mass=me/test 
dir  =2 

[ParticlesN] 

;Particle  emitter  on  dielectric  on  RHS  of  system 
101  shape=Box(xmin-dz,0.0,zmax-2*zbuff-0.25*dz-PMLdepth*dz-Ndielect 
xmax , 1 . 0 , zmax -2* zbuf  f  +0 . 25*dz-PMLdepth*dz-Ndielect ) 
method=BEAM 
temp  =2 . 1/ test 
; temp=10 . 0/test 
energy =0 . 0/test 
106  random=l 

q=-qe/test 
mass=me/ test 
current=2 . 5/ test 
; inject  interval=10 
111  inj  ect  =  l 

f lattop= . 5/f c 
amp2  =0 
tstop=  .  5/f c 
smooth=l 
116  dir=2 

[PlanewaveN ] 

; Generates  Planewave  from  xmin  to  xmax  across  system 
shape=Box (xmin -2*dz , 0 . 0 , zmin+PMLdepth*dz+zbuff , xmax+2*dz , 1 . 0 , zmax 
PMLdepth*dz-zbuff  ) 

121  theta=180 . 0 
phi =0 . 0 
psi =180 . 0 
f requency=f c 

origin=(xmin ,0.0 , zmin+zbuff+dz*PMLdepth) 

126  E0=energy 

[Expert ] 

; cell  mult=12  ; not  sure  what  this  parameter  does  -  not  supported 
in  current  ICEPIC 
r eque  st  ed_max_num_part icles=3e6 
131  self _f i elds = 1  ;ensures  space  charge  is  accounted  for 
[BoundN ] 

; PML  on  LHS  of  system 
method=PML 
depth=PMLdepth 
136  R=1 . 0E-4 
taper =0 
order  =2 
dir  =2 

shape=Box (xmin ,0 , zmin-dz/2 , xmax , 1 . 0 , zmin+dz/2) 
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141 

146 

151 

156 

161 

166 

171 

176 

181 

186 


[BoundN ] 

;  PML  on  RHS  of  system 

method=PML 

depth=PMLdepth 

R= 1 . OE-4 

taper =0 

order  =2 

dir  =2 

shape  =  Box ( xmin ,0 , zmax - dz /2 , xmax , 1 . 0 , ziax  +  dz/2) 

; [ DumpN ] 

; dump_f ormat =GRID 
; dump_plane=l 
; dump_ value =0 
;  nst  op  =0 

[DumpN ] 

dump_f  ormat  =  PART 
dump_interval=DumpInt 
nst  art  =numst  art 
nstop=numstop 
collected_part=l 

[DumpN ] 

dump_f  ormat  =  FIELD 

dump_plane=l 

dump_ value  =0 

fieldflags=l-l-l-2-p-l 

dump_interval=DumpInt 

shape=Box(xmin ,  xmin , (zmax-PMLdepth*dz-Ndielect -2*  zbuf  f - 1*  dz ) -(zmin.  .  . 
+PMLdepth*dz+Ndielect+2*zbuff+l*dz) /2 , xmax , 1 . 0 , zmax-PMLdepth*dz . .  . 
-Ndielect -2*zbuff -l*dz) 
dump _in_ PML =0 

; ; dump_exclude_s t at i c =0  ; not  supported  by  current  version 

dump_name=f ieldSys 
nst  art  =numst  art 
nstop=numstop 

; [ DumpN ] 

; dump_f ormat =CQUADS 

;  shape  =  Box (xmin , xmin , zmin , xmax , 1 . 0 , zmax) 

; nst  op  =0 

[DumpN ] 

dump_f  ormat  =  FIELD 

dump_plane=l 

dump_ value  =0 

fieldflags=l-l-l-2-p-l 

dump_interval=Dump!nt 


80 


191  dump_name=f ieldR 
ns t  art  =numst  art 
nstop=numstop 

shape=box (xmin ,0 , zmin+PMLdepth*dz+(zbuff /2) -0 . 5*  dz , xmax , 1 . 0 , zmin+. 
PMLdepth*dz+(zbuff /2) +0 . 5*dz) 

196 

[DumpN ] 

dump_f  ormat  =  FIELD 
dump_plane=l 
dump_ value  =0 

201  f ieldf lags=l-l-l-2-p-l 
dump_ interval =DumpInt 
dump_name=f ieldB 
ns t  art  =numst  art 
nstop=numstop 

206  shape  =  box (xmin , 0 , zmin+PMLdepth*dz+  (3*  zbuff /2)  -0 . 5*dz , xmax , 1 . 0 , zmin 
+PMLdepth*dz+(3*  zbuf  f/2)+0.5*dz) 

[DumpN ] 

dump_f  ormat  =  FIELD 
211  dump_plane=l 
dump_ value  =0 
fieldflags=l-l-l-2-p-l 
dump_ interval =DumpInt 
dump_name=f ieldT 
216  nst art =numst art 
nstop=numstop 

shape  =  box (xmin , 0 , zmax -PMLdepth*dz - (3*  zbuff /2)  -0 . 5*dz , xmax , 1 . 0 , zmax 
-PMLdepth*dz-(3*  zbuf  f/2)+0.5*dz) 

221 

[DumpN ] 

dump_f  or mat  =  POWER 
dump_plane  =2 

dump_ value  =  zmax -PMLdepth*dz - (3*  zbuff /2) 

226  dump_name=powerT 

shape=Box (xmin , xmin , zmin , xmax ,1.0, zmax) 

nst  art  =numst  art 

nstop=numstop 

231  [DumpN] 

dump_f  or mat  =  POWER 
dump_plane  =2 

dump _ value  =  zmin  +  (3*zbuff /2) +PMLdepth*dz 
dump _ name  =  powerB 

236  shape=Box (xmin , xmin , zmin , xmax ,1.0, zmax) 
nst  art  =numst  art 
nstop=numstop 
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[DumpN ] 

241  dump_f ormat =P0WER 
dump_plane  =2 

dump_value=zmin+ (zbuff /4) +PMLdepth*dz 
dump _ name  =  powerR 

shape=Box (xmin , xmin , zmin , xmax ,1.0, zmax) 

246  nst art =numst art 
nstop=numstop 

[DumpN ] 

dump_f  ormat  = POWER 
251  dump_plane =2 

dump_ value  =  zmax -PMLdepth*dz-2*  zbuf  f -Ndielect -l*dz 
dump _ name  =  power S 

shape=Box (xmin , xmin , zmin , xmax ,1.0, zmax) 
nst  art  =numst  art 
256  nstop=numstop 

[DumpN ] 

dump_f  ormat =PR0BE 
dump_interval=DumpInt 

261  p=((xmax-xmin)/2,0, zmax -2*  zbuf  f -1 . 25*  dz-PMLdepth*dz -Ndielect) 

[DumpN ] 

dump_f  ormat  =  PR0BE 
dump_ interval =DumpInt 

266  p=(  (xmax -  xmin)  /  2 ,0  ,  zmax  -2*zbuff -0. 25* dz-PMLdepth*dz -Ndielect) 

[DumpN ] 

dump_f  ormat  =  J_D0T_E 
avg_ interval =DumpInt 
271  dump_ int erval =1 
nst  art  =numst  art 
nstop=numstop 

shape  =  Box(xmin , xmin , ( zmax -PMLdepth*dz -Ndielect -2*zbuff - 1*  dz ) -(zmin.  .  . 
+PMLdepth*dz+Ndielect+2*zbuff+l*dz)  /2  ,  xmax  ,  1 .0  ,  zmax-PMLdeptli*dz  .  .  . 
-Ndielect -2*zbuff -l*dz) 

276  ; ; [ DumpN] 

; ; dump_f ormat =DENS ITY 
; ; species_num=0 
;  ;  dump.interval =20 
; ; nstart=14*period 
281  ; ; nstop=15*period 
; ; dump_local=0 

;  ; shape  =  Box (xmin , xmin , ( zmax -PMLdepth*dz -Ndielect -2* zbuf  f-l*dz )-(... 
zmin+PMLdepth*dz+Ndielect+2*zbuff+l*dz)/2, xmax, 1.0, zmax - . .  . 
PMLdepth*dz -Ndielect -2* zbuf f -l*dz) 

[DumpN ] 

286  dump_f ormat =RESTART 
dump_ int erval =1 


82 


f  or ce_dump_last =1 
nstart=5*period-2 
nstop=5*period 

291 

[Repeat ] 

Rtype  =  " Dump " 

Rindex=0 

Riaxindex= int ( ( ( zmax - zmin  -  2* PMLdepth * dz ) /6 . 0 ) / dz ) 

296  Rstep=l 

[BeginRepeating] 
dump_f  ormat=VOLTAGE 
dump_ interval = Dump Int 
ns t  art  =numst  art 
301  nstop=numstop 

st art  1  =  ( ( xmax - xmin ) /2 . 0 , 0 . 0 , zmax  -  PMLdepth*dz  -  Ndielect  -  2*... 
zbuf f  ) 

endl =(( xmax - xmin ) /2 . 0 , 0 . 0 ,( zmax  -  zmin  -  2* PMLdepth *dz ) *5 . 0/6 . 0  + 
PMLdepth*dz  -  Ndielect  -  2*zbuff  +  Rindex*dz) 
num_ sample  =4*int((( zmax -zmin-2*PMLdepth*dz)/6.0)/dz)  -  4*(Rindex. 
-1) 

[EndRepeat ] 
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Appendix  C.  Code  for  Calculation  of  Theoretical 
Susceptibility  Curves  for  the  Oblique  Case  in  Matlab® 


%  version  2  incorporates  printing  out  the  yield  based  on  the 
%  phase  -  aver aged  energy,  to  compare  with  the  phase  -  aver aged  yield 
%  oblique  approach  to  multipactor  susceptibility  curves 
%  dimensionless  quantities: 
i  t  =  w  *  t 

7«  v  =  v/sqrt  (Emax/m) 

7«  E  =  q*E/ (m*w*  sqrt  ( Emax /m)  ) 

7«  x  =  w*x/sqrt  (Emax/n) 

7. 

%  this  case  finds  the  phase  averaged  impact  energy,  then 
I  For  these  parameters: 

7. 

7.  f  r eq  =  2 . 45 e  +9  ;  7,  Hz 

7.  w  =  2*pi*freq; 

7,  m  =  9.8e-31;  7.  kg 

i  Emax0=420;  7«  ev 

7.  qe  =  l  .  6  e  -  1 9  ;  7.  eV 

I  Emax  =420  .  *  qe  ;  7o  J 


7.  time  of  1  =  real  time  of  6.496120e-ll  s,  or  0.15915  of 

period 

7. 

7.  velocity  of  1  =  a  real  velocity  of  8.280787e+06  m/s 

7. 

7.  distance  of  1  =  a  real  distance  of  5.379299e-04  meters 

7. 

7«  energy  of  1  =  a  real  energy  of  420  eV 

7. 

7.  efield  of  1  =  a  real  Efield  of  0.78077  MV /m 


7«  therefore,  based  on  Valfells,  et  al  .  figure  5,  we  want  the  ... 

values  (in  the  scale  of  the  figure) 

7.  to  be  as  follows: 

7.  Erf  '  0  to  10 

7,  Edc  ~  0  to  1.8 

7.  which  correspond  to  real  values  as  follows  (figure  units  to  real 
units 

l  conversion  is  2.5:1): 

7.  Erf  ~  0  to  25  MV/m 

7.  Edc  ~  0  to  4.5  MV/m 

7«  which  can  be  translated  to  the  normalized  values  of  (normalized 
to  real  value  conversion  is  1.28e-6:l) 

7.  Erf  ’  0  to  32 

7,  Edc  ~  0  to  5.76 


7. .  .  . 
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%  Debugging 


clear  all ; 
clc  ; 

close  all ; 
seelt  =0 ; 
plotlt  =0 ; 
movieYes=l ; 

u 


%  Global  Parameters  -  these  are  the  only  parameters  that  need  to 
be 

7.  specified 


7. 

eMaxOArg  =420 ; 
ksArg=l ; 
dMaxOArg  =3 ; 
zetaInit=pi/2 ; 
f req=2 . 45e+9 ; 


7.  eV 

7.  the  smoothness  parameter 
7.  secondary  emission 
7.  the  impact  angle 
7.  Hz 


voArg  =  l/ sqrt  (420/2)  ;  7»  initial  velocity 

phiArg  =  pi/2;  7»  emission  angle  (pi/2  is  normal  to  the  surface) 

psiMat  =  deg2rad  (  [5  8  10  20  30  50  75]);  7.  angle  of  obliqueness  (... 

angle  of  Erf  wrt  surface  -  zero  is  parallel) 

ErfMat  =  linspace  (0 , 32 , 50)  7.  the  ranges  of  the  E-fields.  See  ... 

notes  above  for 

EdcMat  =  linspace  (0  ,  5 . 76 , 25)  7.  explanation  of  units 


7.7. 

7. . . . 


7.  From  the  literature  (Valfells,  et  al  -  Plasma  Physics  Jan  98) 
7.  this  is  the  function  way 

7. .  .  . 


7.  phi  is  the  injection  angle  for  the  emitted  particle  (pi/2  is 
normal  to 

7.  the  surface)  Erf 

7.  psi  is  the  angle  of  the  Erf  with  respect  to  surface  (zero  is 
parallel ) 

7.  theta  is  the  phase  angle  of  the  Erf  field 
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ax=@ (Edc ,Erf  ,psi ,t , theta)  -Edc-Erf .  *sin(t  +  theta)  .  *sin(psi)  ; 
ay=@ (Erf , theta, psi ,t)  -Erf .*sin(t+theta) .*cos(psi) ; 


vx=@ (Edc , Erf , phi , psi , t , theta , vo )  vo.*sin(phi)  -  Edc.*t  +  Erf.*sin(... 
psi) . *(cos (t+theta) -cos (theta) ) ; 

vy  =  @ (Erf  , phi  , psi , t , theta ,  vo)  vo.*cos(phi)  +  Erf  *  cos (psi )  .  * ( cos (t+ . .  . 
theta) -cos (theta) ) ; 

x  =  @ (Edc , Erf  , phi , psi  , t , theta , vo )  -Edc. *t  .  “2/2  +  t . * vo . *  sin (phi ) +Erf . .  . 
. *  cos ( theta )  .  .  . 

.  * ( - 1  +  sin(t))  ,*sin(psi)-Erf.*sin(theta)  .*sin(psi)+Erf.*cos(t... 
) .*sin(theta) .*sin(psi) ; 

y=@ ( Erf , phi , psi , t , theta , vo )  t  .  * vo . *  cos (phi )  +  Erf  .  *  cos (psi )  . * ( -t  .  * . . . 

cos (theta) -sin (theta) +sin (t+theta) ) ; 

energy  =  @(Edc ,Erf  ,phi ,psi , t , theta ,vo)  0.5. *(vx(Edc ,Erf  ,phi ,psi  ,t  ,  . .  . 
theta , vo) . ~2  ... 

+vy(Erf ,phi ,psi ,t , theta ,vo) . ~ 2) ; 


7.7. 

7. . . . 


7.  Loop  over  different  values  of  psi 


Erf  Mat  =  Erf  Mat  (2  :  end  ,  1 )  ;  7.  truncate  off  the  first  zero,  since  the... 

solvers 

EdcMat  =  EdcMat  (2  :  end  ,  1 )  ;  7.  won’t  like  zeros  as  the  first  ... 

argument  s 

if  movieYes ==1 

mov=avif ile ( ’ oblique.avi’  ,  ’  f  ps  ’  ,5)  ; 
end 

for  p=l : length (psiMat ) 
psiArg=psiMat (p) ; 


7. .  .  . 


7.  The  middle  loop  starts  here  -  the  goal  is  to  set  EDC  to  a  value 
and  dynamically 

7.  vary  ERF  until  we  find  a  value  of  delta  that  is  close  enough 
to  1  . 

7.  However,  the  two  ’for1  loop  method  will  creep  up  on  it  -  it  . 
will  systematically 

7.  adjust  Erf  and  Edc  and  plot  the  results  with  a  plus  or  minus 
based  on 

7.  growth  or  decay. 
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re  suit  sCell  =  cell  (  length.  (  Erf  Mat  )  *  length  (  EdcMat  ),  3)  ;  °/0  this 

stores  the  Edc  vs  Erf  values  &  line  markers 
amax=length (EdcMat ) ; 
bmax=length (Erf Mat ) ; 

7.  initialize  the  yield  results  matrices 
dPhAvgMat  =0 ; 
dEnAvgMat  =0 ; 

for  a=l : amax 

EdcArg=EdcMat (a , 1) ; 

for  b=l : bmax 

Erf Arg  =  Erf Mat (b  ,  1)  ; 


i  while  abs  (1-deltaAvg)  >tol  7.  one  way  to  do  it  is 
to  dynamically  adjust 

7.  Erf  to  find  the  ... 
value  that  makes  .  . 
delta 

7.  within  a  specified 
tolerance 


7.  Now,  let’s  solve  the  inner  loop  (sampling  50  times 
in  the  phase  of  the  RF 
1  field) 


nmax  =50 ; 

enList  =  zeros  (nmax  ,  1)  ;  7»  initialize  the  storage  ... 

vectors 

tcList=zeros (nmax ,1) ; 
dList=zeros (nmax ,1) ; 
tcGuess  =0.1; 

xtTheta=@(t , theta)  x(EdcArg , ErfArg , phiArg  ,psiArg  ,t  ,  .  . . 
theta , voArg) ; 


7«  a  fast  way  to  generate  titles 
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155 


160  % 

% 

1 

1 


165 


170 


175 


180 


185 


190 


conditions  =  strcat ({ ’ Erf  =  ’ } , num2st r ( Er f Arg )  ,  {  ’  ,  Edc  =... 

’ } , num2str (EdcArg)  ,  ... 

\psi  =  1  }  ,  num2str (rad2deg (psiArg)  )  ,  1 " {\ cir c }  ’  . . . 

); 

for  n=0 : nmax 

thetaArg=2*pi*n/ nmax  ; 
xt=@(t)  xtTheta (t , thetaArg) ; 

•/.  find  the  hoptime,  associated  impact  energy  and  ... 
impact  angle 

tcList (n+1 , 1) =fzero (xt , tcGuess ) ; 

tcGuess  =  tcList  (n  +  1  ,  1)  ;  "/»  update  the  starting... 

value  for  the  fzero  function 
while  tcList (n  +  1  ,  1)  <  0.005 

tcGuess=tcGuess*2; 

tcList (n  +  1  ,  1) =fzero (xt , tcGuess )  ; 

end 

enList (n  +  1 , 1) =energy (EdcArg , ErfArg , phiArg , psiArg  ,  . .  . 

tcList (n  +  1  ,1)  ,thetaArg ,voArg)  ; 
impactAngle  =  abs (vy (ErfArg , phiArg , psiArg , tcList (n.  .  . 
+l,l),thetaArg,  voArg) / . . . 
vx( EdcArg .ErfArg .phiArg .psiArg .tcList (n  +  1 ,1)  , .  .  . 
thetaArg .voArg)) ; 
zet aTi It = at an ( impact Angle ) ; 

•/.  find  El  and  E2  for  this  case  and  the  yield 
[el ,e2 , deltaCase  , eMaxCase]  =  sec_emission(eMaxOArg , .  .  . 
ksArg.dMaxOArg.zetaTilt  ,0) ; 

dList (n+1 , 1) =deltaCase (enList (n+1 , 1) *eMaxCase) ; 

if  seelt  ==1 

disp ([’ iteration  =  1  num2str(n)]) 

disp([’  Hoptime  =  ’  num2str (tcList (n  +  1  ,  1) )] ) 
disp([’  hoptime/period  is  1  num2str (tcList (n  +  1  ,  1)  ..  . 
*2*pi ) ] ) 

disp([’  Energy  =  ’  num2str ( enList (n+1 , 1) )] ) 

disp([’  Phase  (degrees)  =  ’  num2str ( rad2deg ( .  .  . 
thetaArg ) ) ] ) 

disp([’  Yield  =  ’  num2str ( dList (n  +  1  ,  1 ))] ) 
disp([’  El  =  ’  num2str ( e 1 ) ] ) 
disp([’  E2  =  ’  num2str ( e2 ) ] ) 

disp([’  impact  angle  =’  num2str ( rad2deg ( zetaTilt ) ) . .  . 

]  ) 

end 

end 

enAvg  =  mean  ( enList  )  ;  °/„  phase  averaged  energy 

deltaAvg  =  mean  (dList )  ;  7.  phase  averaged  yield 

nrange  =  (0:l: nmax )  ’  ; 
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if  plotIt==l 

scrsz  =  get  (0  ,  ’  ScreenSize  ’  )  ;  7,  get  the  screen  size 

f igure (’ Position  [scrsz (3) /4  scrsz (4) /4  scrsz (3) /2  .. 
scrsz  (4)  /2]  )  ;  70  make  the  figure  in  the  middle  of  . 

the  screen 

plot (nrange , enLi st ) ; 

title ([’ Impact  Energy  vs  Iteration  (2*n*pi/nmax)  for  ’ 
condit ions ] )  ; 
xlabel ( ’ Iteration  ’  )  ; 
ylab el (’Energy’)  ; 

f igure (’ Position [scrsz (3) /4  scrsz (4) /4  scrsz (3) /2  .. 
scrsz  (4)  /2]  )  ;  7«  make  the  figure  in  the  middle  of  . 

the  screen 
plot (nrange , dList) ; 

title ([’ Yield  vs  Iteration  (2*n*pi /nmax )  for  ’  ... 

condit ions ]  )  ; 
xlabel ( ’  Iteration  ’  )  ; 
ylabel ( ’ Yield  ’  )  ; 
end 

[el ,e2,deltaCase , eMaxCase] =sec_emission (eMaxOArg ,ksArg 
, dMaxOArg , pi/2,0) ; 

/  store  the  two  different  yields  in  a  matrix 
dPhAvgMat  =  vert  cat  (  dPhAvgMat  ,  delt  aAvg  )  ;  7,  the  phase  . 

averaged  yield 

dEnAvgMat  =  vert  cat (dEnAvgMat  , deltaCase ( enAvg * eMaxOArg ) ) 
;  7»  the  yield  based  on  the  phase -avg  energy 

7.  display  the  results 
disp([’ERF  =  ’  num2str (Erf Arg) ] ) 
disp([’EDC  =  ’  num2str ( EdcArg ) ] ) 
disp([’Psi  =  ’  num2str (psiArg) ] ) 

disp([’The  phase  averaged  impact  energy  is  1  num2str(. 
enAvg ) ] ) 

disp([’The  phase  averaged  yield  is  1  num2str ( deltaAvg ) 

]) 

disp([’The  yield  based  on  the  phase  -  aver aged  energy  is 
’  num2  str (deltaCase (enAvg*eMaxOArg) ) ] ) 

7.  find  out  if  it’s  growth  or  decay  and  plot  it  with  .. 

the  appropriate 
7.  marker  type 

if  deltaAvg  <  1 

markerType= ’ ob 1 ; 
elseif  deltaAvg  >=  1 
iarkerType= ’+r ’ ; 

end 


7.  store  the  results  in  a  cell  array 
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270 


resultsCell{(a-l)  *bmax+b  ,  l}=EdcArg  ; 
resultsCell{(a-l) *  bmax+b , 2}=ErfArg ; 
resultsCell{(a-l) *  bmax+b ,3}=markerType ; 


end 

end 


*/.  plot  the  results  at  the  very  end 

scrsz  =  get  (0  ,  ’  ScreenSize  ’  )  ;  7.  get  the  screen  size 

f igur e ( ’ Pos it  ion ’  ,  [  s cr sz  (3) /6  scrsz (4) /6  scrsz(3)/1.5  scrsz (4) 
/ 1 . 5 ] )  ;  /  figure  in  the  middle  of  the  screen 

set  (0  ,  ’ defaulttextinterpreter ’  ,  1  latex ’ )  ; 

7.  f  =  figure  ( 1 )  ; 

axi s ( 1 xy  ’  )  ; 

7oSet(gca,  ’Units’,  ’Inches’,  ’  OuterPosition  ’  ,  [0.25  0.25  8  6]) 

i 

'/,axis  (  ’  image  ’  )  ; 

7,  set  (f  ,  ’Units’,  ’Inches’,  ’PaperPosition’,  [0  0  8  6]); 
hold  on; 

for  d  =  l : size ( resultsCell  ,  1) 

plot (resultsCell{d , 1 }  *  ( (le9/freq)*sqrt (400/ eMaxOArg) . .  . 
*(1/1 .28) )  , resultsCell{d,2}*((le9/freq)*sqrt (400/.  .  . 
eMaxOArg) * (1/1 . 28) ) , . . . 
r esult sCell {d  ,  3} )  ; 

end 

xlabel ({  ’  ’;’  ’;’$E_{dc,  actual}  [MV/m]  \times  (f/1  GHz)~{\... 

mbox{-}  1} (E_{max  ,  0}/400  eV) "{\mbox{-}  1/2}$  ’  }  ,  ’ FontSize ’ . . 
,12)  ; 

ylabel ({ ’ $E_ {rf ,  actual}  [MV/m] \times  (f/1  GHz ) “ [\mbox { -} 1} ( E_ 
{max  ,0}/400  eV)"{\mbox{-}l/2}$’;’  ’;’  ’},  ’FontSize’, 12); 

title ([’ Susceptibility  Curve  for  $\psi  =$  ’  num2str ( rad2deg ( . .  . 
psiArg))  ’$“\circ$  \&  $\delta_{max , 0}  =$’  num2str ( dMaxOArg) 
]  ,  .  .  . 

’ FontSize ’  ,  14)  ; 

[hx,hy]=format_ticks(gca,  ’$$’  ,  ’$$’)  ; 
hold  of  f  ; 

savef igs ( [ ’ oblique_PSI = ’  num2str ( r ad2deg ( ps i Arg ) )  ’deg’]); 

if  movieYes==l 

winsize  =  get(gcf  ,  ’Position’)  ; 
winsize ( 1 : 2)  =  [0  0]; 

F(p)=getframe(gcf .winsize) ; 
mov=addframe(mov,F(p)) ; 

end 

close  all ; 

f igur e ( ’ Pos it  ion  ’  ,  [  s cr sz (3) /4  scrsz (4) /4  scrsz (3) /2  scrsz  (4)  . . 
/  2  ]  )  ;  7.  figure  in  the  middle  of  the  screen 
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hold  on; 

plot(dPhAvgMat  ,  ’  b  ’  ,  ’LineWidth’  ,2)  ; 

275  plot ( dEnAvgMat  ,  ’  r  ’  ,  ’LineWidth’  ,2)  ; 

title ([’ Yield  for  $\psi  =$  ’  num2str (rad2deg (psiArg) ) 
circ$  \&  $\delta_{max ,0}  =$  ’  num2str ( dMaxOArg ) 

’ FontSize ’  ,  14)  ; 

xlabel({’  ’  ;  ’  ’j’Edc/Erf  pair’},1 FontSize’, 14); 
ylabel ({’ Yield FontSize 14)  ; 

280  [hx,hy]=format_ticks(gca,  ’$$’  ,  ’$$’)  ; 

legend (’ $\theta$  averaged  $\delta$ ’ , ’ $\delta$  based  on  $\... 
theta$  avg  energy’); 

savef igs ( [ ’ y i eld_PSI = ’  num2st r ( r ad2deg ( ps i Arg ) )  ’deg’]); 

close  all ; 

end 

285 

if  movieYes==l 
mov=close (mov) ; 
end 
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Appendix  D.  Input  File  for  the  Oblique  Case  in  ICEPIC 


Listing  D.l.  Sample  Input  File  With  ip  =  5° 

[Def  ault  s ] 

2  ENERGY  =  1 . 5  e6 
EDC  =0  .  1 25  e6 

s inps i =0 . 0872  ;  sin(psi)  is  a  constant  for  each  simulation 
;  this  is  based  on  fichtl  ’s  original  input  files  and  is  meant  to 
expand  it  to  apply  to  oblique 
;  fields  (8  Dec  08) 

7 

[Variables] 
c0=299792458 . 0 
pi =3 . 141592653589793 
eps0=8 . 85418781762E-12 
12  me  =9 . 1093897E-31 
qe  =  1 . 60217738E-19 
test  =  1 
f c=2 . 45E9 
Iamda0=c0/ f c 
17  ez=EDC 

energy = ENERGY 
dz=lamda0/3000 

PMLdepth=int ( ( . 025* lamdaO ) / dz  +  .  5) 

Nx  =  50 
22  Nz  =300 
xmin=0 
xmax=Nx*dz 
zmin=0 
zmax  =  Nz  *  dz 

27  dt  =0 . 99  *  dz  /  (  cO  *  sqrt  (2)  ) 
period=int (1/fc/dt) 
max step  =  15*  period 
dumpint  =  period/30 
numst  art  =  0* period 
32  numstop=3*period 

hoptime=(2*me/qe/ez)*sqrt (2*qe*420/me) *0.8509035245 
hopdist=sqrt (2* qe *420/ me) *hoptime 
Ndielect=10*dz 
print  "Hop  time=",  hoptime 
37  print  "Hop  distance=",  hopdist 

print  " TimeSt eps /Hopt ime = " ,hoptime/dt 
print  "Hop  distance /dz=" ,  hopdist/dz 

[Cartesian] 

42  ig_dir  =  (0 ,1,0) 

def ault_e = (0 . ,0. ,-ez) 

; f lattop= . 5/f c 
; amp2  =0 
t  st  art  =0.0 
47  tstop=10*period 
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cO  =  l 

cl =( ENERGY /- EDC ) *sinpsi 
f requency=f c 
ramp=5*dt 
52 

[XGridl] 

range =Uni form (xmin -3*dz , xmax+3*dz , Nx+6) 

[YGridl] 

57  range  =  Unif orm (0 , 1  ,  1) 

[ZGridl] 

range  =  Unif  orm ( zmin , zmax , Nz ) 

62  [Symmetry] 

; symmetry  =  (xmin,-l.  ,zmin-10*dz, xmax ,2.  , zmax+10*dz) 
periodicl=(X , xmin , xmax ) 

;  sym  mode  =  (  1 . 0  ,  0 . 0  ,  0 . 0  ) 


67 

[Time] 
dt  =  dt 

step_max=maxstep 
courant_value= . 95 
72  kill_if _below_part s =5 
kill_if_above_parts=2e6 

[ShapeN ] 

shape=Box(xmin , xmin , zmin ,xmax ,xmax ,zmax) 

77 

[ShapeN ] 

; Dielectric  on  RHS  of  system 

shape=Box(xmin ,0.0 , zmax -4*dz-PMLdepth*dz-Ndielect ,xmax ,1.0, zmax -4* 
dz-PMLdepth*dz) 
material=dielectric 

82 

[ShapeN ] 

; Dielectric  on  LHS  of  system 

shape=Box(xmin ,0.0 , zmin+4*dz+PMLdepth*dz ,xmax , 1 .0 , zmin+PMLdepth*dz 
+Ndielect+4*dz) 
material=dielectric 
87 

[MaterialN] 
name=dielectric 
epsilon=l . 001 

92  [ParticlesN] 

; Secondary  emitter  on  dielectric  on  RHS  of  system 

shape  =  Box (xmin-dz ,0.0,zmax-4.25*dz-PMLdepth*dz-Ndielect  ,xmax  ,1.0,  . 

zmax -3 . 75*dz-PMLdepth*dz-Ndielect ) 
method = SECONDARY 
sec_coefficient=3 
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97  sec_energy=2 . 1/ test 
sec_threshold=l/  test 
sec.refl =0 . 03 
sec.scat =0 . 07 

sec_energy_max_yield=420/ test 
102  sec_ks=l 

q=-qe/test 
mass  =me /test 
dir  =2 

107  [ParticlesN] 

;Particle  emitter  on  dielectric  on  RHS  of  system 

shape=Box(xmin-dz,0.0,zmax-4.25*dz-PMLdepth*dz-Ndielect  , xmax  ,  1 . 0  , 
zmax  -  3 . 75*dz-PMLdepth*dz-Ndielect  ) 
method=BEAM 
temp  =2 . 1/ test 
112  energy=0 . 0/test 
random=l 
q=-qe/test 
mass  =me /test 
current =1000/ test 
117  in j e ct _ int erval = 100 
f lattop= . 5/f c 
amp2  =0 
tstop=  .  5/f c 
smooth=l 
122  dir  =2 

[PlanewaveN ] 

;  Generates  Planewave  from  xmin  to  xmax  across  system 
shape=Box ( xmin -2*  dz , 0 . 0 , zmin+PMLdepth*dz+2*dz , xmax+2*dz , 1 . 0 , zmax- 
PMLdepth*dz-2*dz) 

127  theta=180.0 
phi =0 . 0 
psi =180 . 0 
f requency=f c 

origin=(xmin ,0.0 , zmin+2*dz+dz*PMLdepth) 

132  E0=energy 

[Expert ] 

; cell_mult =12  ; not  sure  what  this  parameter  is... 

r eque  st  ed_max_num_part icles=3e6 
137  self _f i elds =0 

[BoundN ] 

; PML  on  LHS  of  system 
method=PML 
142  depth=PMLdepth 
R= 1 . 0E-4 
taper =0 
order  =2 
dir  =2 


94 


147 

152 

157 

162 

167 
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192 


shape=Box (xmin ,0 , zmin-dz/2 , xmax , 1 . 0 , zmin+dz/2) 


[BoundN ] 

;  PML  on  RHS  of  system 

method=PML 

depth=PMLdepth 

R=1 . 0E-4 

taper =0 

order  =2 

dir  =2 

shape  =  Box  (xmin  ,0  ,  zmax  -  dz  /2  ,  xmax  ,  1 . 0  ,  zmax  +  dz/2) 


[DumpN ] 

dump_f  ormat  =  FIELD 
dump_plane  =  1 
dump_ value =0 
fieldflags=l-l-l-2-p-l 
dump_interval=dumpint 

shape=Box(xmin , xmin , (zmax-PMLdepth*dz-Ndielect -5*  dz ) -(zmin+. 
PMLdepth*dz+Ndielect+5*dz)/2,xmax , xmax ,zmax-PMLdepth*dz-. 
Ndielect -l*dz) 
dump _in_ PML =0 

; dump_exclude_s t at i c =0  ; not  supported  by  current  version 

ns t  art  =numst  art 

nstop=numstop 

; [ DumpN ] 

; dump_f ormat =CQUADS 

;  shape  =  Box  ( xmin  ,  xmin  ,  zmin  ,  xmax  ,  xmax  ,  zmax  ) 

;  nst  op  =0 

; [ DumpN ] 

; dump_f ormat =P0WER 
; dump_plane  =2 

; dump _ value  =  zmax -PMLdepth*dz-3*dz 
; dump.name =powerT 

;  shape  =  Box  ( xmin  ,  xmin  ,  zmin  ,  xmax  ,  xmax  ,  zmax  ) 

; nstart=22*period 
; nstop=25*period 

; [ DumpN ] 

; dump_f ormat =P0WER 
; dump_plane  =2 

; dump_ value = zmin +3* dz+PMLdepth*dz 
; dump.name =powerB 

;  shape  =  Box  (xmin  ,  xmin  ,  zmin  ,  xmax  ,  xmax  ,  zmax) 

; nstart=22*period 
; nstop=25*period 

; [ DumpN ] 

; dump_f ormat =P0WER 


197  ; dump_plane =2 

;  dump_value  =  zmin  +  l*dz  +  PMLdepth.*dz 
; dump_name =powerR 

;  shape  =  Box  (xmin  ,  xmin  ,  zmin  ,  xmax  ,  xmax  ,  zmax) 
; nstart=22*period 
202  ; nstop=25*period 
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[ DumpN ] 

dump_f  ormat  =  PROBE 
dump_ int erval =1 

p= ( ( xmax -xmin ) /2 ,0 , zmax -5 . 25*dz-PMLdepth*dz-Ndielect ) 
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[ DumpN ] 

dump_f  ormat=PR0BE 
dump_ int erval =1 

p= ( ( xmax -xmin ) /2 ,0 , zmax -4 . 25*dz-PMLdepth*dz-Ndielect ) 
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[ DumpN ] 

dump_f  ormat  =  J_D0T_E 
avg_interval=1000 
dump_ int erval =100 
nst art =22* period 
nstop=25*period 

shape=Box (xmin , xmin , (zmax-PMLdepth*dz-Ndielect -5*  dz ) -(zmin+ 
PMLdepth*dz+Ndielect  ;+5*dz)/2,xmax  ,xmax  , zmax-PMLdepth*dz- 
Ndielect -5*  dz ) 


222  ; [DumpN] 

; dump_f ormat =DENSITY 
; species_num=0 
; dump_ int erval =20 
; nstart=14*period 
227  ; nstop=15*period 
; dump_local =0 

; shape=Box (xmin , xmin , (zmax-PMLdepth*dz-Ndielect -5*  dz ) -(zmin+ 
PMLdepth*dz+Ndielect+5*dz)/2,xmax  ,  xmax  ,zmax-PMLdepth*dz-. 
Ndielect -5*  dz ) 
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Appendix  E.  Mathematica®  Notebook  for  the  Approximation 
of  the  Theoretical  Power  Deposited  on  the  Surface  of  a 

Dielectric 


CODE  BEGINS  ON  NEXT  PAGE 
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Verboncoeur  Inputs  for  ICEPIC  results 


■  This  notebook  accumulates  the  results  presented  by  Verboncoeur 
in  IEEE  transactions  on  Plasma  Science,  Vol.  28,  No.  3,  June  2000 
pp  529-536.  We  want  to  develop  a  working  form  for  interpreting 
the  results  of  the  ICEPIC  runs. 

We  assume  that  our  inputs  will  include: 

Basic: 

N  =  number  of  pseudo  particles 
L  =  electrode  separation 
At  =  integration  time  step 
p  =  particle  weighting 

Case  specific  inputs 

ERf  =  Rf  field  amplitude 
f  =  frequency  of  Rf  field 
T  =  temperature  of  emitted  electrons 
or  alternatively 
vt  =  thermal  velocity 

The  final  input  is  the  ejection  current  density,  J0 


J0  =  VtT  n0e  vt/2 

Since  we  have  assigned  vt,  the  sole  parameter  is  n0.  Our  purpose 
is  to  lay  out  the  methodology  to  make  this  assignment  and 
estimate  the  power  deposited  on  the  surface  of  a  dielectric  in  the 
space  charge  case.  Therefore,  we  must  first  start  by  finding  the 
value  of  n0by  integrating  Verboncoeur's  equation  (15). 


In[45]:=  Clear  [  "Global"  *"  ] 
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ln[46]:= 


ln[57]:= 

ln[61]:= 

Out[61]= 


dx  =  3 . 05910671429  *  10  A (-5) ;  (*the  length  in  the  x  direction*) 
dy  =  l*A-2 ; 

zmax  =  0.0016;  (*the  total  length  in  z  for  the  interaction 

space  of  our  system*) 

vt  =  9 . 92*A5 ; 
qe  =  1 . 602*A-19; 
eps  =  8.8542*A-12; 
me  =  9.11*A-31; 

«p[n_]  =  n*qeA2/  (eps*  me)  ; 

nmacro = 1*A5;  (*from  the  ICEPIC  run  -  ice.stat  file  *) 
ratio = 1 . 2384*A3;  (*this  is  the  ratio  of  macroparticles  to 

electrons  -  fixed  for  a  given  inject 
current  *) 

ntot  =  nmacro  *  ratio;  (*this  gives  us  the  total  number  of  real 

electrons  in  the  system*) 

nzideal  =  ntot  /  (dx*dy);  (*  this  is  the  ideal  total  number  of 

real  electrons  in  particles/mA3  - 
we  need  to  set  the  integral  of 
Verboncoeur ' s  equation  15  equal 
to  this  *) 

Print ["Ideal  Nz  =  ",  nzideal] 

Ideal  Nz  =  4 . 04824  xlO14 

int2  =  Integrate [1 /  (1  +  a  * z)  A 2 ,  {z,  0,  zmax}]; 

al  =  Vn*qeA2/  ( eps  *  me  *  vt  A  2 )  ; 
f  [n_]  :  =  Evaluate  [  (int2  *  n  /  .  a  -*  al)  ]  ; 
f [nn] ; 

f  [5 . 29  *  10 A 20]  (*this  is  the  number  that  should  match  up  with 

nzideal*) 

4 . 043  x  1014 
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We  consider  the  y-component  of  the  velocity  and  assume  the 
injected  value  of  this  component,  vOy,  is  zero. 


f  =  2 . 45  10  A9;  Ef  rf  =  3 . 0  x  10  A  6;  n  =  5 . 29  *  10  A  20;  shield  =  1; 

(*  Shield  ==  1,  shieling  case  for  G,  Shield  ==  0  , 
no  shielding  *) 

(*  the  n  value  here  is  nO  from  the  above  calculations*) 
w  =  2  7r  f ; 

qe  =  1 . 6  x  10  A  (-19)  ; 
me  =  9.1xlOA  (-31)  ; 
eps  =  8.8542  x  10A  (-12)  ; 
toev  =  1.  /  (1.6xlOA  (-19))  ; 
dx  =  3. 0591067142  9*  10  A  (-5)  ; 
dy  =  1; 

cO  =  2 . 99792  *  10  A  8; 

vrf  =  qe  *  Efrf  /  («  me)  ;  kT  =  2 . 8  /  toev;  vt  =  V 2  kT  /  me  ; 

Print ["The  Rf  electric  field  velocity  is  vrf  =  " 

,  vrf,  "  [m/sec] "] ; 

Print["The  thermal  velocity  is  vt  =  ",  vt,  "  [m/sec]"]; 

(*  Plasma  frequency  *) 

wp  [n_]  :  =  -\J n  qe  A  2  /  (eps  me) 


(*  Efx  E  field  due  to  surface  charge  *) 

Efx[n_]  :=  2  kT  n  /  eps  ; 

Print ["x  E-field  ",  Efx[n]] 

Print  ["x  E-field  velocity  ",  qeEfx[n]  /  (mew)] 

(*  cycle  time  or  bounce  time  Eqn[19]  *) 

z  [v_,  n_]  :=  V 2  7r  Exp[vA2  /  (2  vtA2)  ]  Erf  [Vv  A2  /  (2  vt  A2)  ]  j 
wp[n]  ; 


(*  Shielding  function  *) 
Gshield[U_,  shield_]  :=  shield*  (l 


NIntegrate 


[cos  V 2 


7r  U  Exp  [yA2]  Erf 


[y]  ] 


Exp [ -  2  y  A  2 ] 
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{y,  0 ,  Infinity}  J  j  +  (1  -  shield)  *  (1 .  -  Exp  [  -U  A  2  ]  )  ; 

(*  Current  density  *) 

jZero  [n_]  :  =  V^r"  n  qe  vt  /  2  .  ; 

Print [ "Current  Density  jZero  =  ",  jZero [n],  "  A/m2"] 

Print ["Ratio  of  Plasma  Frequency  to  Rf  frequency  ",  wp[n]  /  a]; 

(*  Average  Impact  energy  *) 

ai[U  ]  :  =  me  *  vrf  A  2  /  2  *  Gshield  [U  ,  shield]; 

Print  [  "Average  Impact  energy  in  ",  ai  [w  /  wp  [n]  ]  *  toev,  "  eV  "] 

Print  [  "Maximum  impact  energy  =  ",  2  me  *  vrf  A  2  toev,  "  eV"  ]  ; 

(*  Average  Power  deposited  *) 

powerDeposited [n_]  :  =  ai  [w  /  wp[n]  ]  *  jZero  [n]  /  qe; 

Print ["Power  deposited  ",  powerDeposited [n] ,  "  W/mA2"] 

Print [ "Average  Power  Deposited  ",  powerDeposited [n]  * 
dx*dy,  "  W"] 

uenergy  =  0 . 5  *  eps  *  Ef  rf  A  2 ; 

poynting  =  cO  *  uenergy; 

ptheory  =  dx  *  dy  *  poynting; 

Print [ "Percentage  of  power  deposited  ",  powerDeposited [n]  *  dx 

The  Rf  electric  field  velocity  is  vrf  =  3.42652  xl07  [m/sec] 

The  thermal  velocity  is  vt  =  992278.  [m/sec] 

x  E-field  7 . 31656  x  106 

x  E-field  velocity  8.3568xl07 

Current  Density  jZero  =  7.4431 xlO7  A/m2 

Ratio  of  Plasma  Frequency  to  Rf  frequency  84.2183 

Average  Impact  energy  in  3.08009  eV 

Maximum  impact  energy  =  13355.5  eV 

Power  deposited  2.29254xl08  W/mA2 

Average  Power  Deposited  7013.14  W 

Percentage  of  power  deposited  1.91927% 
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